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Preface 


This book is, more or less, the script of a lecture series entitled “Phase-Field Theory 
and Application,” which I give at Ruhr University Bochum (RUB) in my department 
Scale-bridging Thermodynamic and Kinetic Simulation (STKS). It strongly rests on 
my personal experience during many years in academic and applied research with a 
special focus on "phase field." The idea is to combine a very brief presentation of the 
materials-related mechanisms alongside the theoretical background of the phase- 
field theory and, last but not least, to introduce numerical solutions. The lecture 
series is designed for graduate students and newcomers to phase-field theory. There 
are many aspects of phase field that will only be covered briefly in this script. At 
my university, we are continuing this lecture series with an advanced course in the 
format of a discussion seminar. We hope that we can continue with this in a second 
volume of the book. 

The 12-lecture course, as held at RUB, is condensed to seven lectures in this 
book. An eighth lecture, *Quantum Phase Field;" is added here that I do only in the 
advanced course: Just free your mind and enjoy! 

Exercises are recommended at the end of each lecture, and suggestions for further 
reading are given. Examples are added where appropriate. A short introduction to 
the phase-field code OpenPhase academic [1] is added as a "tutorial" with two 
examples as a beginning to your own research. 


Bochum, Germany Ingo Steinbach 
June 2022 
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Part I 
The Phase-Field Method 


Chapter 1 A 
Introduction FECA 


1.1 What Is a Phase Field? 


A phase field, in the context of these lectures, is a scalar field in space and time 
that indicates the local state of matter as characterized by its phase state. This may 
not help much for now: do not worry, we will proceed step by step, explaining all 
expressions. Firstly, the present section will explain what a phase field is—i.e., what 
do we understand by this expression? Secondly, in Sect. 1.2, we will investigate 
which types of materials problems can be examined using phase fields. In Sect. 1.3, 
the historical background of phase field will be detailed; perhaps you will come 
back to this part at the end of the lectures, at which time you should be better able 
to reference the different historical aspects of phase field comparing to their current 
state of development. In the last section of this introductory chapter, Sect. 1.4, we 
will focus in particular on the scale of the materials problem: the mesoscopic scale. 
Let us begin! 

The “phase” of matter, a lump of atoms within a reference volume, denotes the 
state of crystalline order between the atoms. It is a central concept in physics and 
a central element of the phase-field theory. The order of atoms—i.e., their position 
in a crystalline lattice: FCC, BCC, or other crystalline structures— defines the phase 
state. In addition, the amorphous, liquid, and gaseous states—i.e., the absence of 
crystalline order—are characterized as phases. If the phase state is a gas, we may 
be referring to an atmosphere outside of a solid-state sample or a pore inside the 
sample. Liquids are mainly considered as a melted phase in connection with a 
solid crystal. This can be water and ice, or molten iron in connection with already- 
solidified iron and slag. Further types of phase are characterized by magnetic or 
ferroic order, superconductivity, or plasma. The phase field, a field in space and time, 
indicates the state of a materials point (in space and time), whether solid, liquid, gas, 
plasma ferromagnetic, etc, It tells nothing about other properties of matter at this 
point (in space and time), like its temperature, pressure, the composition of elements 
or molecules. Of course this information is also needed. But it is not enough to 
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characterize the material at this point (...). We need to know it's phase state, which 
is enclosed in the local value of the “phase field." 

In the physics literature, the phase state is characterized by a so-called "order 
parameter,” which is normalized between 0 and 1. We will use the letter $ to 
represent this parameter. The term “order” relates to “crystalline order," which is 
different between different phases. For example, a crystalline solid will be stable at 
low temperatures, while its melt will be stable above the melting temperature. An 
important element of the concept of phases is that there is a discontinuity in their 
order; the order does not change gradually if, for example, the temperature changes: 
there is an abrupt change in the atomic order (at least for phase transformations of 
the first kind, which we will be considering here). This discontinuity is also reflected 
in discontinuities in the properties of the phases, such as their elasticity, viscosity, 
thermal or electrical conductivity. In other words: the phase state can be determined 
uniquely from the material properties of the piece of matter we are investigating. 

Another important aspect is that the phase state needs not to be stable. We also 
consider metastable or unstable phases. This means that two materials with the same 
temperature, pressure and composition may be in different phase states. In other 
words, to characterize the state of a material not only pressure, temperature and 
composition is needed, but also information about the atomic order: its phases state. 

To make this concept more transparent to materials scientists, let us discuss a 
metallic alloy, say a binary aluminium-silicon alloy, regarding possible phases at 
different temperatures and compositions. Figure 1.1b shows the phase diagram of 
this alloy. This is called a “eutectic” phase diagram, as it has three stable phases 
in different regimes of alloy composition c and temperature T, as indicated in the 
figure: liquid, FCC aluminum, and FCC silicon (diamond lattice) on the far right- 
hand side of the diagram is almost pure silicon, which is not shown. These regimes 
in which a phase is stable are also termed “phase fields" in alloy thermodynamics. 
Don't be confused. They are fields in temperature, composition, and pressure, not 
in space and time as we will use the term in the context of these lectures. But both 
usages of "phase field", of course, are somewhat related. 

In-between these stable regions of individual phases, we have the so-called “two- 
phase" regions, where none of the adjacent phases can be stable. The reason for this 
can be read from the Gibbs energies of the system as function of composition [see 
Fig. 1.1a]. For a given temperature and pressure, the Gibbs energy functions Gy 
and Gg of the two phases, o and £, are displayed schematically. The minimum of 
the total Gibbs energy is determined by the so-called double-tangent construction, 
where the common tangent touches the curves at the phase compositions c4 and cg. 
With the fraction of phases fy = 1 — fg, we define the total Gibbs energy Gtotal: 


Gtotal = faGa + fgGg. (1.1) 


The total Gibbs energy is minimal along the tangent line between cy and cg, 
weighted with the fractions fy and fg. Each individual phase with composition 
€ (cy < c < cg) will have a higher energy. Therefore, the material with a 
nominal composition co in the two-phase region must decompose in two phases 
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Ca Cp Csilicon 
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Fig. 1.1 (a) Schematic of the Gibbs energy curves of a binary alloy in two different phase states, 
a and f, for a given temperature To. The common tangent determines: (i) the phase compositions 
Cq and cg of an alloy with nominal composition within the two-phase region, and (ii) the minimum 
Gibbs energy of a two-phase mixture with this nominal composition. (b) Linearized Al-Si phase 
diagram, indicating the stable composition of cy and cg for the given temperature To 


so that the material will reach a lower state in the Gibbs total energy. We say 
that the composition region between cy and cg is "forbidden" for both phases. 
There is an "energy barrier" in the Gibbs energy of the material that separates it 
into different phases. Other properties of the different phases—such as elasticity, 
thermal conductivity, and density—also show a similar behavior, and they are 
clearly different in both phases. 

All of this is well known to most of our readers, but we repeat it here to 
perform the following thought experiment. The upper panel of Fig. 1.2 displays a 
cross section through a two-phase mixture between liquid (white phase) and FCC 
aluminum (black phase) close to equilibrium. A line scan of the local composition 
is performed, and the corresponding composition c(x) is displayed in the lower 
panel of this figure. This switches from cy to cg and back again. It is now an easy 
exercise to normalize this composition to determine the phase fractions as a function 
of space, as also displayed along the line scan: 


SE. fg —— io fgml- fe (1.2) 


f a= , 
Cy — Cf CB — Ca 

The fractions fy and fg are usually used to characterize a phase mixture as a 
whole; the fractions are not considered as fields in space and time. Therefore, we 
use a different notation: we use a(x), the “phase field,’ which varies in space 
between 0 and 1. Note that we will also consider interfaces between phases as 
“diffuse,” i.e., with a finite width. Therefore, the graph in Fig. 1.2 is not step-wise: 
it is “smeared out,’ as will be discussed in detail in Sect. 1.3.1. The normalized 
function "phase field" is also called an "indicator function" because it indicates 
where a phase is present in space. In the same way as it is defined from the varying 
phase concentrations in (1.2), it could be defined from the elasticity or the density of 
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Fig. 1.2 Scheme ofa 

solid-liquid phase mixture 

(solid in black, liquid in 

white) close to the 

equilibrium of the 

compositions within the 

individual phases. The scan 

evaluates the local 

composition, which varies 

between the solid 

concentration cy and the 

liquid concentration cg. On Øsotid 
the left-hand axis, the value 

of the phase field @solia IS l 
displayed, with 8 indexing 

the solid phase 


the individual phases. We can also map the phase field back to the local composition, 
elasticity, and density if we know the local phase-field value and the properties of 
the individual phases. 

The phase field is a field in space that indicates the position at which we find a 
special phase a, Øy (x) = 1, or do not find this phase, øy (x) = 0. Then we let this 
evolve in time (!): 


fa > Pax) > by (x, t). 


We can see that a multi-phase material can be characterized regarding its local 
phase state. We can use the standard concepts of alloy thermodynamics to determine 
the local phase fractions in a small reference volume containing, say, a thousand 
atoms to be large enough to characterize the phase state reliably. Then, we will find 
either the values of 1 or 0, or if our reference volume is intersected by an interface, 
the value will be between 1 and 0. Interfaces may move, and phases can change. 

The content of the following chapters involves establishing how to determine the 
evolution of phases in complex materials under various conditions. 


1.2 What Is the Purpose of Phase Field? 


Generally speaking, the phase-field method uses a set of partial differential equa- 
tions to describe a material's problem with evolving microstructures. Let us first 
elaborate on what we mean by “evolving microstructures.” The microstructure of a 
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multicrystalline material—and we will mainly speak about the microstructures of 
metallic materials—is the distribution and shape of different grains in that material. 
These individual grains have many attributes: their orientation with respect to a 
reference orientation; their composition, pressure, and temperature; along with 
important crystalline defects such as dislocations, stacking faults, and vacancies. 
If two neighboring crystalline grains are of the same phase but have different 
orientations, they form a “grain boundary.” Consequently, they are considered as 
different elements of the microstructure and will be denoted by different indicators, 
that is, by different phase fields. Two grains of different phases form a "phase 
boundary,” and a pore within a solid grain forms a “surface.” 

We will describe all these different cases using phase fields. This means that a 
phase field marks a region of space in which the material can be i) characterized 
uniquely by its phase state and orientation or ii) as an interface (where the phase 
state is undetermined). Practically speaking, the phase field $ (X), a field in three- 
dimensional (3D) space x, has, by convention, the value $4 (X) = 1 for the material 
point at position x belonging to the phase o indexed by $4. For I > $4(X) > 0 
it is an interface, for $4 (X) = 0 it is another phase than a. Other conventions can 
be found in the literature, but this is the prevailing one. If there are only two grains 
in our material, e.g., a small inclusion in a matrix, we will skip the subscript œ and 
index the inclusion by $(x) = 1 and the matrix by $ (x) = 0 (or vice versa). The 
phase boundary is marked by an intermediate value 0 — $(x) < 1. Treating the 
field as continuous, or “diffuse,” as we do in phase-field theory, the boundary has to 
have a finite width n. 

Such a construction is very helpful for characterizing a polycrystal, a multi- 
phase material, or any general static microstructure. We want, however, to allow the 
microstructure to change in time, $ (x) — $ (x, t). Now we have to distinguish two 
cases as follows. (i) The microstructure is just “deforming,” i.e., we find a mapping 
of space coordinates X(fj) — X(t2). In continuum mechanics, this is called the 
mapping from the “reference frame" at time f, to the “actual frame" at time t2. In 
the second case, (ii) the microstructure is "evolving"; it may be growing, swelling, 
or shrinking. It may have been born in the time t; > 1? ("nucleated,” as we say in 
materials science), or it may have died, i.e., it has fully transformed to a different 
phase. Thus, it cannot just be referenced in two different frames. In general, one 
material point in an evolving structure transforms to a material point of a different 
kind (e.g., liquid to solid). No physical movement of matter has to be involved in this 
process. We will speak in these chapters about case (ii), evolving microstructures. 
We will in addition allow them to deform (see Chap. 7), but this is a side remark. 

The example shown at the end of this chapter is a dendrite growing in two 
dimensions (2D). In this example, which we undertake as an exercise, only the 
interface between the phases (solid and liquid) “moves.” This means that a material 
point, which had been liquid in the starting condition, changes its phase state to 
solid. However, no material is moving, neglecting the effects of shrinkage during 
the phase transformation and motion of the solid and the melt. For small changes 
in the shape of the dendrite, we could also treat this problem using finite elements. 
Then, the elements would have to be swelling or shrinking in both phases. However, 
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in the general case, if we start from full liquid and simulate to full solid, we will 
need so-called adaptive finite elements, which continuously adapt to the shape of 
the dendrite. This is not impossible: it has been demonstrated by several authors, 
particularly in 2D simulations, but also in 3D (e.g., by Provatas et al. [16] and 
Schmidt [17]). However, this is numerically expensive, particularly in 3D. The 
phase-field method offers an elegant and numerically efficient way to treat this 
problem of evolving, growing, or shrinking microstructures. Aside from using 
adaptive finite elements, alternative approaches include the cellular automata and 
level set methods (see Further Reading). 

Having described a materials problem—e.g., solidification or coarsening of a 
microstructure during heat treatment—as a moving-boundary problem using a set 
of partial differential equations, we will also need external boundary conditions and 
initial conditions. Furthermore, we will need a good database for physical input data 
and constitutive models for all the properties of the bulk phases and boundaries. 
Last but not least, we will need good solvers and sound numerics. All of this will be 
touched on during these lectures. It should be noted that the scientist applying the 
method will have to be patient: a good numerical solution to a problem may take 
weeks or months. . .. 


13 History of the Phase-Field Method 


1.3.1 Microscopic Phase Field 


The phase-field method can be split into two branches with very different histories, 
interpretations, and intentions. The first branch is often referred to as “order- 
parameter theory," the “microscopic phase-field method,” or the “time-dependent 
Ginzburg—Landau theory.” This rests on grounds in thermodynamics with a special 
focus regarding interfaces. Commonly, van der Waals is noted as being the father 
of this branch; he rationalized from general considerations that a diffuse interface 
between two phases will be more likely than a sharp interface [23], although 
at this time, the existence of the atomistic structure of matter had not yet been 
established. As a next step, Landau introduced the concept of an order parameter 
into the thermodynamic description of materials. Ginzburg must be referenced 
as introducing gradients of the order parameter into the concept, representing 
interfaces or phase boundaries. Other stepping-stones in the development are the 
Cahn-Hilliard theory of spinodal decomposition [3] and Khachaturyan’s theory of 
microelasticity [8]. Wheeler, Boettinger, and McFadden introduced a first model 
of alloy solidification in 1992 [24], and this combines a Cahn—Hilliard model of a 
material with a miscibility gap and a phase-field model of the second type, which we 
will call a “mesoscopic phase-field model” below. A compilation of both branches 
with notable contributions, by far incomplete, is shown in Fig. 1.3. Before giving 
the historical details of the second branch—mesoscopic models—let us provide an 
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Phase-Field Theory ~ 1994 
Korteweg Plapp 
~ 1895 2005 e 
Langer Kim (KKS) 
- 1978 1999 
(J Fix Steinbach 
~ 1983 1996 
Caginalp Karma 
@ ~ 1986 - 1994 oe 
[2 Mesoscopic e wm» Fracture 
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Fig. 1.3 Two branches of phase-field theory—microscopic and mesoscopic—highlighting 
important steps in their development. Today, the branches converge to a common understanding 
accepting the strengths of each approach. We stop this history about 20 years ago, so newer 
developments in phase-field theory are not included, such as models of fracture [6, 14, 18] or 
quantum phase fields [11, 21] 


outline of the problem that urges us to extend microscopic models, along with its 
possible solution. 


1.3.2 The Problem of Scale 


All phase-field models—and this must be very clear—agree in terms of their general 
structure; they start from a thermodynamic functional density f (6$, Vo, T,c,€,...) 
with three terms of different character: 


f- (V9 + L9! 0 «mij. T, NEA (1.3) 


Here, we use the notation of Kobayashi [9]. The first and second terms relate to 
interface or grain-boundary contributions. This is easy to see, since the first term 
vanishes in the bulk when ¢ = const, i.e., in the bulk of the grains ¢ = 1 ord = 0 
(in the convention we shall use throughout this script). For these conditions, the 
second term—the so-called “double-well potential”—also vanishes by construction. 
Both terms represent a positive energy penalty within the interface, 0 < $ < 1, 
which is related to the interface energy (see Chap. 2). The last term in (1.3) relates 
to the bulk energy difference between phases. This is a function of temperature T, 
composition c, strain e, and other material states such as magnetism. Charges may 
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be added, but those cases are not treated in this lecture series, and they are little 
investigated in phase-field theory. 
We can also write down the functional in the form: 


f =Z [Por e nea - 67] + aec. T. ese...) (1.4) 


where o is the interfacial energy, 7 is the interfacial width, and Ag is the energy 
difference between phases, the function m in Eq.(1.3) above. Now, if all the 
parameters of the equations—e, y, m, o, n and Ag—are constants, we can find 
a direct one-to-one correlation between the two representations, as will be shown in 
Chap. 2. 

We can take for granted that the relation between m and AG works, since it 
is not phase field specific: it is standard thermodynamic Gibbs energy difference 
between different phases. The relations between e, y, o, and 7 have the following 
problem: if € and y are determined from an underlying microscopic theory, as 
well as the interface energy o, the interface width 7 is uniquely fixed. If all 
parameters are correct, we will find, as in real materials, that the interface width 
n % lnm. This is good if we are investigating materials at the nanometer scale 
(as scientists do with microscopic phase-field models); it is not so good if we 
want to investigate microstructures at the micrometer scale, because it will not be 
feasible to treat the diffuse interface of a phase field with a width of 1 nm in a 3D 
simulation of micrometer-sized grain structures. How can this dilemma be solved? 
We seek a theory that can be matched to a so-called "sharp interface" theory (see 
Chap. 4). In this theory, the interface width has no physical importance. To match 
a diffuse-interface method (the phase-field method) to the theory of such a sharp 
interface, we need to remove the influence of the interface width on all quantities 
with physical meaning. We seek a theory that is agnostic regarding the interface 
width of a real material: the interface width shall be scalable for convenience of 
numerical simulation. To say this clearly: it is frustrating to throw away the physical 
insights of a phase field to predict the structure of an interface; it is, however, the 
great success of mesoscopic phase-field models to make quantitative predictions of 
microstructural processes at larger scales possible! 


1.4 Mesoscopic Phase-Field Model 


The mesoscopic phase-field model! stems from a very different route than thermo- 
dynamics. It is rooted in wave mechanics, in particular “traveling-wave solutions" 
called “solitons” in the physics literature. The investigation of these phenomena 
dates back to Korteweg, a Dutch mathematician in the nineteenth century. It is 


l The phrase “phase field” was adopted by Jim Langer in 1978 [13] for a special phase 
transformation: solidification. 
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(a) (b) 


Fig. 1.4 Schematic representation of solitonian waves. The left-hand panel displays a wave pulse 
Ø(E) of width 7 along the space coordinate x traveling with speed v; this remains self-similar in 
the moving coordinate & = x — vt. The right-hand panel shows a wave front; this is an integral 
form of the wave pulse. The width 7 and moving coordinate & are the same as for the wave pulse 


reported that he observed a single wave packet, caused by a boat clutching against a 
channel wall in Amsterdam, his hometown. The wave travelled for a long distance 
without changing its shape. He and his collaborators wrote down a non-linear wave 
equation for such a phenomenon in the shallow waters of rivers or ocean shores. 
This is known today as the Korteweg-de Vries (KDV) equation; it is of third order 
in the space derivative, and its solution is displayed in Fig. 1.4a. 

The special feature of this solution is self-similarity in its local coordinate 
$ = x — vt with the velocity v of the wave packet, withstanding moderate 
outer perturbations. This property, as realized by Langer, can be used to propagate 
interfaces by solitonian waves in a numerical simulation: the shape of the wave 
front is not part of the solution, i.e., it can be eliminated from the solution. The wave 
packet, Fig. 1.4a, hereby does not distinguish between either side of the packet, right 
or left. Therefore, in phase-field models, the integral form of the KDV equation 
is used as a second-order equation with an asymmetric solution, as displayed in 
Fig. 1.4b. This is more or less exactly what we call today phase field: the minimum 
solution of (1.3) or (1.4). 

In the physics literature, both the symmetric and the anti-symmetric solutions 
are termed “solitons”: self-similar traveling-wave solutions of non-linear wave 
equations, in contrast to periodic solutions of linear wave equations. There are 
further aspects of this connection between linear and non-linear wave equations that 
relate to the fundamental understanding of matter (see Chap. 8). In particular, soliton 
solutions can be used for quantization if applied to an appropriate wave function, 
as already pointed out by Olsen in 1974 [15]. This "quantization of phase field" 
has been applied to studying scale formation in the universe, [22]. Enthusiastically 
speaking, quantization of phase field opens the road to a new understanding of the 
physical world: space, time, and matter are unified in a wave-mechanical description 
that is consistent with both quantum mechanics and thermodynamics (see further 
reading, “Solitons and quantum phase field,” on this topic). 


12 1 Introduction 


Returning to the classical phase field in materials science, the idea is to combine 
the equation for the solitonian front, which travels with a given velocity v, with a 
transport equation in the bulk phases ahead or and behind the wave. The transport 
in the bulk phases ahead and behind, of course, is influenced on the one hand by the 
wave itself. On the other hand, the state of the phases ahead and behind will affect 
the wave velocity wave v. For the classical soliton, the front velocity is a constant 
such that the traveling wave is a plane wave; however, now in our applications, the 
velocity is a vector in 3D space, varying at each position of the wave front dependent 
on the interaction of the wave front with its surroundings: v > v(x, t).2 The 
equation used to couple both phenomena is the empirical Gibbs- Thomson equation. 
This relates the velocity of the front linearly to the kinetic undercooling of a front 
with mobility M^ and the front normal n: 


ù =ñ M? (c*k + Ag). (1.5) 


The Gibbs-Thomson equation considers capillarity with the interface stiffness 
c* and the local curvature of the front «; i.e., the kinetic undercooling is expressed 
as the curvature-corrected Gibbs energy difference between the bulk phases Ag. The 
Gibbs-Thomson equation is not dependent on the interface width 7: it represents 
a sharp-interface model at the mesoscopic scale, where the interface can be 
approximated as a discontinuity between the phases. The Gibbs-Thomson equation 
(1.5) can be replaced by an appropriate phase-field equation (see Chap. 2). Since the 
width of the phase-field function in the transition between phases has no analogue 
in the Gibbs-Thomson equation, it can be chosen (within certain bounds) for 
numerical convenience. This will be treated in detail in Chap. 4. 

It was also realized early on [1, 2] that the curvature « is considered inherently 
in a phase-field description of a curved interface. This will be discussed in detail in 
Chap. 3. 


1.4.1 Applications 


The phase-field method now has many applications, from biophysics to astro- 
physics. Its main field of application is materials science, and metallurgy in 
particular. As described above, we distinguish the two main approaches: micro- 
scopic models, which are applied mostly to solid-state transformations (the school 
of Armen Khachaturyan); and mesoscopic models, which were first applied to 
solidification (the school of Jim Langer). Today, there are well-established methods 
for combining both branches for specific materials problems; these range, in 


? Throughout these chapters, we will only use vector notation where it is needed for understanding. 


1.4 Mesoscopic Phase-Field Model 


the context of metallurgy, from solidification and heat treatment to in-service 
degradation until failure. A new branch has arisen with phase-field models of 
fracture [7, 14, 18]. This will, however, not be covered in this lecture series. One 
common aspect of all these applications is the evolution of the microstructure. 
Evolving structures define a special class of moving-boundary problems, in which 


phase-field models offer elegant and efficient solutions. 


Example—Dendritic Growth 

Dendritic solidification occurs if a solid crystal grows from an undercooled 
melt. As a result, the interface becomes unstable, forming tips and branches. 
It then develops into a “tree-like” structure: hence the term “dendrite.” 


Equiaxed Dendrite 

This example is related to Kobayashi's dendrite in 2D [9, 10], and it is an 
exercise that we undertake in the class. One very good result from a former 
student is displayed in Fig.1.5. The “100-line code" (see Appendix A.1), 
written in C++, evolves the phase-field variable and temperature field with 
adiabatic boundary conditions. One circular grain is placed at the bottom of 
the 2D simulation zone at the start of the simulation, and the system size is 
300 x 600 grid points. The phase field $ is displayed for different time steps. 


Further Reading 


Long-range elastic, electrostatic, and magnetic interactions connected to 
domain structure evolution during structural, ferroelectric, and ferromag- 
netic phase transformations, a recent review: [4]. 

Fracture: [6, 14, 18]. 

Level set: [5, 20]. 

Hydrodynamic interpretation of Phase Field: [12] 

Solitons and quantum phase field: [11, 19, 21]. 
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(a) (b) 


(c) (d) 


Fig. 1.5 Dendritic solidification of a pure substance into an undercooled melt. (a) t = 0.05 s. (b) 
t = 0.25s. (c) t = 0.75s. (d) t = 1.5s 
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Chapter 2 A 
Analytics om 


2.1 The Problem of Propagating a Wave Front on a 
Numerical Grid 


In this chapter, we will recapture the analytical background of traveling-wave 
solutions in 1D with two of the most commonly used separating potentials: the so- 
called “double-well” and “double-obstacle” potentials. This solution, also known 
as the antisymmetric soliton solution (see previous Chap. 1), is the basic reason 
that phase fields work so well in a numerical solution: the solution is quite robust 
against perturbations, which are inevitable in numerical simulations. Of course this 
resembles the physical fact, that the wave front is robust against any perturbation— 
wind, a bird etc.—in reality. The front self-stabilizes while traveling over the 
numerical grid, just like a solitonian wave on the surface of water. 

Let us start with an example. We define a wave front $ (x, fo) at time fo in one 


dimension x: 
3 (1am (5°) 
O(x, to) = = | I — tanh | — : (2.1) 
2 n 


We know the wave front from Fig. 1.4b: the solitonian wave, but now defined 
between 0 and 1, as we do nowadays in phase field. In the following we take this 
convention that the left-hand side of the wave is $ = | and its right-hand side @ = 0. 
The front is centered around x = 0, and the width 7 is a measure of the distance 
between $& = 0.05 and $ = 0.95. This "cutoff" is necessary since the hyperbolic 
tangent converges to O or 1 only at infinity, and we need a finite measure of the 
interface width. 
We now want to propagate this wave with constant speed vg. One writes: 


d(x, t) = —— = —v. (2.2) 
X X 
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You may try to solve this equation numerically in 1D (as we do in the class), 
applying the profile (2.1) as a starting condition. However, this will simply not work 
for propagating the wave for a reasonable distance while maintaining its profile! 
Equation (2.2) is called a hyperbolic transport equation, and it is prone to large 
numerical instability up to shock-waves, and we do not want this at all. The equation 
thus has to be regularized using an appropriate method. We will not go into the 
details of numerical regularization approaches in hydro-dynamical wave mechanics, 
but we will modify the equation for better solvability. We want to solve (2.2) exactly. 
The only things we can do to the equation are: (i) add 0, or (ii) multiply by 1. We 
will do the first. We add: 


8? 72 1 
TAE. ELT o (5 $). Q3) 


^ 9x2 


Does this look strange? Let us do the derivation. We will prove that (2.3) has the 
solution (2.1) and that it is self-similar in time for a propagating planar wave, a wave 
in 1D. Therefore, it adds 0 to the original Eq. (2.2). We add it: 


(x,t) 204 29 in 
Ox 


8? 72 1 8 
E 79 DIE 2E an (2.4) 


EFT? ax 


This is the famous “phase-field equation” we will be dealing with in the rest of the 
lectures.! Adding this strange 0 Eqs. (2.3)-(2.1) changes the type of the equation 
from hyperbolic to parabolic. It becomes a “diffusion equation,” which is easy to 
solve in a discretized manner on a numerical grid. Let us go through this step by 
step. 


2.2 Equation of Motion for the Phase Field 


We start with the phase-field functional density, Eq.(1.3). A functional, by its 
mathematical definition, is a mapping of an arbitrary set of functions (functional 
densities), within a given definition domain, onto the real numbers. These functional 
densities, functions in space and time, are very difficult to compare. If we map 
functions onto a scalar, the comparison becomes easy. The functions can then be 
ordered on a linear coordinate indicating whether they are “larger or smaller.” For a 
human, sadly speaking, “money” is an often-used measure; this measures the value 
of all goods that can be bought. In metallurgy, we like to use the free energy of the 


! The equation is given here without physical prefactors, which will be used to multiply the 0 later. 
Why? To make you curious for Chap. 3: because in 3D, the “0” gets a physical meaning! 
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System because it is a thermodynamic functional of the state variables temperature 
and pressure, two quantities that are easily controllable in lab experiments. Other 
functionals are useful under different conditions. In thermodynamics, we use 
Legendre transformation to relate these functionals uniquely. Nevertheless, different 
functionals should be used for different applications. 

To introduce “time,” we start from an entropy functional S. The second law of 
thermodynamics demands that entropy increases with time: 


dS 0485 


(2.5) 


Here, the symbol ô denotes the “variational derivative" of a functional, i.e., of a 
function of functions. This will be elaborated below in some detail. For a rigorous 
derivation, see Kiselev, Shnir, and Tregubovich's book [2]. For now, let us simply 
take the message that Eq. (2.5) can easily be fulfilled if 


ó óF 
i x 3 or pu X : (2.6) 
ot óQ ot bd 


This relation is termed the “Clausius—Duhem relation" in continuum mechanics. 
It is not a “physical law,” but it is the simplest possible Ansatz to ensure the 
positivity of entropy production; it works well and is generally accepted. It is 
termed the “relaxation Ansatz" because it has only a first derivative in time and 
no accelerations. Since the entropy S enters the Gibbs free energy F with a “—” 
sign, we accept that: 


F = Jo d?xf is the total free energy within the 3D domain Q with the free- 
energy density f. We have to work on the functional as an integral over the domain 
of interest since the free energy density f of a phase-field model is a function of 
gradients in $, V, which are, by construction, not defined on a special point in 
space but need at least two different points. Let us treat the three parts of the free 
energy in (1.3) separately: 


fi = 3009, (2.7) 

h= rga - 9y., (2.8) 
39?—29? 

få = m@,T,c,€,...)=h(@)Ag(T,c, €,..) = — ^. C,6,...). (2.9) 


We start with the easiest one: f2. We are not afraid of the phase-field function 
(x,t); we treat it as a normal variable at the point in space x and time t 
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under consideration, and it should not depend on the values of the field at other 
locations regarding the function f». This will come later. Therefore, we can treat the 
variational derivative as a normal partial derivative. We calculate:? 


6 F5 EL A 
Te. yọ — Q) (5-9): (2.10) 


As mentioned before, this part is only active in the interface. It changes its sign at 


d= 1, 1.e., in the center of the interface. It is the first derivative of the so-called 
double-well potential, also called the Landau potential, Qu — oy}. Therefore, we 
call it the potential contribution. 

In the third term, f3 (2.9), we have already separated the phase-field-dependent 


part h($) = ue the so-called coupling function, from the physical part, 
i.e., the Gibbs-energy-density difference between the phases as a function of 
temperature, composition, strain, and other parameters, Ag(T, c, e, ...). Here, we 
will take this as a constant Ag = Ago. Because this part also does not depend on 
non-local contributions, i.e., gradients of the phase field, the functional derivative 
can again be treated as a normal partial derivative: 


6F3 9f 
bb. 


= o(1 — $)Aqgo. (2.11) 


As we can see, the special form of (2.9) has been chosen such that the free-energy 
difference acts only within the interface when 0 < $ < 1. Ago is called the “driving 
force,” since it makes a positive or negative contribution to the evolution of the phase 
field &, driving it to grow or to shrink. 

The first part (2.7) is more involved because it contains a gradient in $, (Vo)’, in 
1D, EA We will perform the variational differentiation in 3D so as not to increase 
the confusion of the problem with the dimensionality of the variational derivative. 
The 1D case works identically. Applying the chain rule, we find: 


DOE ( f drevo) = 5 | Px ever» Q.12) 
w 38 Ma" 2 ~ 86 Jo l 


The difficulty lies in finding a way to treat the variation 6(V@) of the gradient 
operator. It helps to realize that all derivatives—regardless of whether they are 
functional derivatives 0$, partial derivatives dø, or partial derivatives in space 


? Note that the common convention in physics and mathematics literature about variational 
derivatives ó violates the otherwise accepted convention that a symbol å should be dimensionless. 
F is an extensive (absolute) quantity, while f is an intensive quantity (defined per volume). The 
variational derivative removes the volume integral and the volume increment dx. 
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Vø—are linear operators and are therefore commutable. What we will do is 
commute the operators ó and V by a partial integration in space: 


f d?x (eV$8(V$)) = — f d?x (ev?) 5 + boundary integral, (2.13) 
Q Q 


óF| d 
39 = —eEV*Ø. (2.14) 


The boundary integral is omitted mostly with the reasoning that the interface 
should not touch the boundary of the domain. This will not be the case in most 
applications, so boundary conditions should always be treated with care. This topic, 
however, will only be touched on in some numerical exercises. The most important 
factor of the derivation is the change of sign in the gradient contribution (2.12) due 
to the partial integration. In the free-energy functional, both interface contributions 
(2.7) and (2.8) are positive penalties against forming interfaces; they will counteract 
with opposite signs in the equation of motion: The Laplace- or diffusion-operator 
will act as to smear out the profile, while the potential operator will act as to sharpen 
it, as we will see in detail. Collecting all terms defines the equation of motion of the 
phase field, the so-called phase-field Eq. (2.4), and it now has physical constants in 
Kobayashi's notation. It is derived from the Clausius-Duhem relation (2.6) with the 
proportionality constant t, which is called the relaxation constant. 


ð 1 
ro? Leve yeu DE 4) HI — $)Ago. (2.15) 


ot 


For completeness, a common formal replacement of the functional derivative 
for theories with gradient contributions (quantum mechanics in general) has the 
following form, known as the “Euler-Lagrange equation,” in which V $ is treated as 
a symbolic entity in the denominator of the first differential operator: 


S r-| * xl 
sp aV 239] 


2.3 Traveling-Wave Solution for the Double- Well Potential 


A “traveling-wave solution," as introduced in Chap. 1 (Fig. 1.4), has the special 
feature of self-similarity in time (for constant driving force Ago and velocity vg); 
i.e., we can define a coordinate & = x — vot to describe this solution in its own frame 
moving with constant velocity with regard to a resting frame. We derive this solution 
as the minimum solution of the functionals (1.3) or (1.4). Of course, the condition 
of self similarity ae = 0 and the minimum energy condition are related, because 
a minimum solution requires the vanishing of the first derivative of the phase-field 
ô? x SF — Oin Eq. (2.15). Here, we will not try to derive the solution by applying 


ot ^ 59 
some scheme to solve partial differential equations, but we will prove that the given 
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solution from the literature does the job!? The solution has already been introduced 
in (2.1) for the initial time step t = 0. We change x > & = x — vot to find the 


traveling solution: 
1 3(x — 
$6.05 (: am ( - 2). (2.16) 
n 


To prove that this is a solution of (2.15), we have to compute the first and second 
derivatives in the space coordinate x: 


a 6 
å Se ad), (2.17) 
x n 
(09) 72 1 
Je a PG $). (2.18) 


We find: 


2 1 12 1 
eve = yo eque mc (Ze = ») $(1— $) (5 = 4) 
n 2 
72 
=0 for "d Ly. (2.19) 


This is the first important result, and all phase-field models agree: the minimum 
solution of $ for a phase-field functional (1.3) with Ag = 0, or with a self-similar 
solution and arbitrary but constant velocity v, demands an interface width 


772€ 
n= |—. (2.20) 
Y 


Inserting this relation back into the phase-field Eq. (2.15), only the driving force part 
remains, and we have: 


ap n ð$ 
— =-¢(1 — d)Ago = = — Ago- 2.21 
GET $(1 — $)Ag0 = 27; Ago (2.21) 
We compare this with the Gibbs-Thomson equation for a moving planar 
interface, i.e., where the capillarity term is 0, and the interface mobility is M^: 


vo = M? Ago, (2.22) 
og op 

— = M? Ago. 2.23 
at 3, BH (2.23) 


5 A reference to the first scientist who derived this solution is difficult. It is not a particular solution 
for a phase field; it is very general. 
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This relates the relaxation constant r to the interface mobility M?: 


7] 


CSS 
6M? 


(2.24) 


Finally, we have to relate € and y to the interface energy o. The interface energy 
(in units FA must come out if we integrate the free energy density (in units FA 
in the normal direction through the interface: 


+00 
o=f a[o + tea - 9] 


oo 


-Foo 18 y 
[. ax (Bet ra oy 


I 2 2 
3 $^(1— 9) 


ll 
qe. 
g + 

8 

A 

Se 


1 
= f d Y Pa- 9» (2.25) 
l 1 
=+ dé ny é0-9 (2.26) 
1 € 
- a = o (2.27) 
og 


We have used the first derivative of @ with respect to x, 57, (2.17) twice, along 
with the relation between 7, y, and 7 (2.20), and we have substituted the integration 
over x by the integration over $ in (2.25). To summarize this derivation, we insert 
the relation between the model parameters in Kobayashi's notation e, y, and t and 
the physical parameters M^, ø, and 7 into (2.15) to arrive at the final phase-field 
equation with a double-well potential: 


ap 
ot 


I 2, 72 1 6 
= M" j0 | Vo md PA "a $)^go(. | Q28) 


Going back to our original problem, to solve the propagation of a wave front 
with constant speed vo (2.2), we now can prove (2.3) just by inserting the second 
derivative (2.18). In fact, both terms cancel if the front is in the right contour. The 
first derivative (2.17) then motivates the ansatz for the driving force m(@) « Žo x 
$ (1 — $) (2.9). Now we see that the hyperbolic transport Eq. (2.2) is turned into a 
parabolic equation with the second derivative in the space coordinate x. This type 
of equation—also termed a diffusion equation—is well behaved and easily solved 
on a computer. You can try this out using the program with which you tried to solve 
the hyperbolic Eq. (2.2). 
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2.4 Interpretation of the Phase-Field Equation 


Let us now have a closer look at the phase-field Eq. (2.28) to determine the effect 
of each term on the phase-field profile. The first term is the Laplacian, or diffusion 
operator. This smoothens any profile, ensuring a smooth phase-field profile between 
0 and 1. For the right-hand hyperbolic tangent profile, it contributes negative 
increments above & > I and positive increments below & < 1 as indicated in 
Fig. 2.1a. The second term stems from the double-well potential. Its first derivative 
changes its sign at @ = 1 because the potential function is symmetric around 5. Now 
we see that the potential contribution has negative increments where the Laplacian 
has positive increments, and vice versa [Fig. 2.1b]. 

This competition between the two contributions guarantees the stable phase- 
field profile! If a numerical solution deviates from the correct profile, the 
contributions will not balance but push the contour back to the correct contour. 
Later, in Chap. 3, we will see that this balance is also violated for a curved interface 
in more than one dimension. This will be used to consider capillarity. Furthermore, 
Ag will not be a constant in real applications: it will be positive or negative, and 
it will vary from point to point in space and time. It will depend on transport of 
temperature, solute, and momentum, as will be explained in later chapters. 


(b) 


L — $ 


-f' pot i 


Fig. 2.1 The competition between spreading and contraction of the phase-field contour. Descrip- 
tion see text 
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2.5 Phase-Field Equation and Traveling- Wave Solution for 
the Double-Obstacle Potential 


As a general comment, we shall keep in mind that the potential function (2.8), which 
is called the double-well potential, is not a unique choice for any system. It was first 
introduced by Landau to explain the behavior of a ferromagnet at the critical point, 
which is called the Curie point in this system. Only at this critical point, which is a 
second-order phase transformation, is it justified to truncate the potential describing 
the interface energies—the so-called Landau potential—to this simple form (for 
details, see [6]). In phase field, however, we are dealing with a first-order phase 
transformation such as solidification or solid-state phase transformations, which are 
far away from a critical point. We also stick to the notion of a mesoscopic phase-field 
method (Sect. 1.4) whereby phenomena inside the interface cannot, and need not, 
be interpreted physically. Therefore, we have some freedom to choose a potential 
function in a different form than (2.8). 

For reasons that will be explained in Chap. 6, the multi-phase-field method, as 
implemented in OpenPhase [3], applies the so-called double-obstacle potential: 


bo y20 
SNØ (2.29) 


The "obstacle" is realized by the absolute signs, which flip the negative branches 
of the function $ (1 — $) to positive values.* Figure 2.2 shows a comparison between 
the double-well and double-obstacle potentials. 


Fig. 2.2 The “double-well” f 
and “double-obstacle” j 
potentials 


4 Note that these absolute signs are omitted in several publications for better readability, but they 
are definitely necessary for the model to work. In the numerical calculations, they are realized by 
a sharp cutoff against negative values of $ < 0 as well as against values larger than 1,4 > I. 
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As we can see, the graphs of the two potentials are similar. They are determined 
such that: (i) the free energy goes to infinity f — -Eoo for $ — +00; (ii) they have 
two minima at f = 0 for $ = 0 and $ = 1; and (iii) the region of the so-called 
potential barrier between 0 and 1 has an area that is i the interface energy o. Other 
potentials can be used that fulfill above criteria, but they are little discussed in the 
phase-field literature. As long as we are speaking about a two-phase system, they 
should not show another minimum. More complex potentials, so-called “multi-well 
potentials" are used for multi-phase systems (e.g., [1, 4]). 

The minimum solution of the free energy with the double-obstacle potential 
(2.29) is: 


for x— vot < —2, 


= 5 sin (2520) for —ł<x-—vt< +}, (2.30) 


pO, t) = I< 


O n= = 


for x — vot > +4. 


As before, we derive in the transition region 0 < $ < 1 (do this as an exercise!): 


0 
-— T V6 - 9. (2.31) 
x n 

(38)? 7? I 

rc GO (2.32) 


The final task is to find a proper coupling function mP corresponding to the 


first derivative in $, (2.31), of the phase-field profile (2.30) for the double-obstacle 
potential. We only give here the function and leave the proof as an exercise: 


“oo = 1 [es — D/$(1— $) + Sedap — »| ^ (2.33) 
- 2 80- . 


Repeating the analysis for the double-well potential, we find the relations 


between the model parameters €P°, yP°, and «P9 and the physical parameters n, 
c, and M?: 
cbO zx? ePO Do 72 n 
Zu o = —— = i M? ==. 2.34 
" y Po 8n ” 8 8 PO Tem 


Collecting all these pieces, we define the free-energy density in physical notation 
for the double-obstacle potential: 


8 
fae DE +77 |o(1— DI + Ag. (2.35) 
qm] 
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The normalization of the interface contribution 2 is necessary to fulfill the 


condition that the integral of this term in the normal direction through the interface 
becomes the interface energy o [cf. (2.26)]. We also chose the notation where o 
is divided by a length, the interface width 7, to emphasize that this contribution is 
an energy density like Ag. The term in the brackets is now a dimensionless contour 
function indicating the interface position. In the physical parametrization, the phase- 
field equation with the double-obstacle potential reads: 


3 * f] 
s = Mt lav - (5 2] “Vou Daso}. (2.36) 


2.6 Gibbs-Thomson Limit of the Phase-Field Equation 


Summarizing all contributions, we can return to the underlying physical equation 
of interest: the Gibbs- Thomson equation (1.5). We have already done this in 3D, 
relating the velocity of the interface v to the rate change of the phase field with the 


normal vector to the interface n = a 
OP OX; = s 
— — = Vov, 
d Ox; ot d 
v= Lu (2.37) 
IV 


The curvature of the interface x, which is given by the strange expressions in 
Fig. 2.3, will be derived in Chap. 3. 


Double-obstacle Double-well 
— n? (RGA n) : sr end m|. 
= TT n Ivo| n " Lm 6T 7 [Vo] + a 
v= M, [o K+ gol v = M, [0" k + Ago) 


Fig. 2.3 The relations between the phase-field equations in Kobayashi's notation and the Gibbs— 
Thomson equation in physical notation for the double-well and double-obstacle potentials 
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2.7 Exercises 


Exercise 
Prove the first and second derivatives of the traveling-wave solutions for the 
double-well and double-obstacle potentials (2.17), (2.18) and (2.31), (2.32). 


Exercise 

Prove that the coupling function in (2.33) does the job, i.e., that the derivative 
of this function with respect to ¢ is proportional to the first derivative of $ in 
the case where a double-obstacle potential is used. 


Exercise 
Prove the relation (2.34) between the model parameters and the physical 
parameters for the double-obstacle potential. 


Exercise 
Check the relations between the model parameters as indicated in Fig. 2.3. 


References 29 


Further Reading 
* Appendix A in the review “Phase-field models in materials science" [5], 
where another potential function, the “top hat" function, is introduced. 
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Chapter 3 A 
Capillarity PEEN 


3.1 Curvature of a Phase-Field Contour 


In the previous chapter, the meanings of individual terms, the gradient and potential 
contributions of the free-energy functional [Eqs. (2.7) and (2.8), respectively], and 
their functional derivatives were discussed. In the free-energy functional, both terms 
contribute 5 of the interface energy o. In the phase-field equation of motion, either 
for the double-well potential (2.28) or double-obstacle potential (2.36), they cancel 
to 0 in ID for the correct analytical solution. If the numerical solution deviates from 
the analytical solution, either because one has started from wrong initial conditions 
or because of numerical errors, the two terms act to push the interface to the right 
solution, i.e., they stabilize the interface contour. For a curved interface in 3D, it has 
already been mentioned that these terms do not cancel! They have the property to 
evaluate “curvature.” This will be proven in this third chapter. 

To summarize: the phase-field equation has two important properties: (i) to 
propagate the interface with a velocity related to the bulk free-energy difference Ag 
between different phases such that the energetically lower phase will grow and the 
other phase will shrink; (ii) to correct this transformation process for “capillarity” 
related to the local curvature of the interface. The latter also considers interface- 
energy anisotropy, which is a dominant effect in many metallurgical processes. This 
will be detailed in the second part of this chapter. Let us first prove the relation 
between the surface terms in the phase-field equation and curvature. 

The formal mathematical definition of the curvature « of any vector field is the 
divergence of the normal vector 7i: 


(3.1) 
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Of course, we evaluate this for our phase field ¢. Applying the divergence to the 


numerator and denominator of [ea gives: 
Vo 


k= E Vo Vg ve| 

We| [Vg] 

VÆR. ve Veg Vo] (3.2) 
Vo Ive 

= 7 [vg -i V|Vo]] 

= F [v un (3.3) 

eee Ive 2 e Ji (3.4) 
Vo n" 2 


We have inserted the definition of 7 into Eq. (3.2). This projects the divergence 
V Vo onto the normal through the interface, so this operator can be replaced by 


the second derivative in the normal direction in Eq. (3.3). We know this second 
derivative from the 1D analysis, Eqs. (2.18) or (2.32) for the double-well or 
double-obstacle potentials, respectively. In the last Eq. (3.4), the relation for the 
double-obstacle potential is inserted. The relation for the double-well potential 
works analogously. We see that this is exactly the expression used in the previous 
chapter for the evaluation of the capillarity term in the Gibbs-Thomson equation 
(cf. Fig. 2.3). 


3.2 Anisotropy of Interface Energy 


In the previous chapter, dealing with the analytic relation of a phase-field equation 
and the phenomenological Gibbs-Thomson equation for interface motion including 
capillarity, all entities—i.e., the phase-field mobility M?, the interface energy c, 
the interface width 5, and the thermodynamic driving force Ag— were taken as 
constants. This is, in general, not the case; sometimes it is a strong oversimpli- 
fication. The mobility and the interface energy are anisotropic functions of the 
orientation of the interface with respect to the orientation of the adjacent phases. 
One distinguishes so-called "inclinations" of the interface with respect to the 
crystallographic orientation of individual crystal phases, and the misorientation 
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of these phases if they are crystalline solids. In general, we have to specify 
three misorientation angles and two inclinations: five angles in total. This poses 
a challenge for determining these functions, but it also generates challenges for 
phase-field models in terms of coping with these functions. 

The mobility M? is relatively easy; it is simply a scalar proportionality factor 
between time f and the so-called kinetic driving force acting on the interface, i.e., 
the "capillarity-corrected bulk free-energy difference.” It will, in general, vary from 
point to point in the interface, because the misorientation and inclination angles of 
a curved interface vary along the interface. Depending on the physical model used 
for the mobility, it may cause numerical problems regarding time-stepping schemes 
(because it relates to time). However, from a conceptual point of view, it is easy 
and does not need further discussion here. In contrast, the interface energy is more 
difficult. 


3.2.1 Interface-Energy Anisotropy: Phenomenological Picture 


Interface-energy anisotropy is a conceptual challenge. For phenomenology, we will 
only discuss the simple case of a crystalline inclusion in an amorphous matrix, say a 
liquid. Then, the dependence of the interface energy on misorientation (three angles) 
vanishes. We will only consider two dimensions; then, only one inclination angle 
0 with respect to the crystal axis 09 remains (see Fig. 3.1). The Gibbs-Thomson 
equation [see Chap. 1, Eq. (1.5)] becomes: 


v — iM (o*k + ^g). (3.5) 


Here, o* = o + o”, with o" = 250, is called the "interface stiffness,” i.e., the 
resistivity of the interface to bending. o” is also called the “Herring torque,” as it 
was introduced by Herring in 1952 [3]. This tells us that the interfacial energy of a 
precipitate can be reduced by two separate mechanisms: (1) that the interface area 
is reduced; and (ii) that the interface is turned into a low-energy inclination. In fact, 
for the equilibrium shape of the precipitation surrounded by its melt, the interplay 
of both mechanisms has significant consequences. Equilibrium here is defined as 
the chemical potential over the system being constant. Furthermore, the precipitate 
neither grows nor shrinks, i.e., v = 0. Since the precipitate has a curved interface, 
there is a potential jump of o*« between the bulk energy of the precipitate and the 
matrix, the so-called “capillarity pressure." 
The Gibbs-Thomson equation (1.5) becomes (v — 0): 


0= ok + Ag 


o*x = —Ag = constant. (3.6) 
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Fig. 3.1 (a) Sketch of the equilibrium shape of a crystal with strong interface anisotropy. The 
preferred crystallographic axis is rotated by 6o from the x axis. The surface normal at one specific 
point at the rounded corner is tilted by 0 from the preferred axis. (b) The interface energy o and 
interface stiffness o* as functions of 0 — 4 


If o*« = constant all along the interface, then o* has to be inversely proportional 
to x, or, because 0* is a material constant, « has to be inversely proportional to o*. 
Where the stiffness is high, the curvature of the interface has to be low. If the 
stiffness is low, the interface can be easily bent, and we have edges or corners. 
We also know that the interface energy is high at edges and corners, so the interface 
area must be small. As a take-away message: where the interface energy is high, 
the stiffness is low, and vice versa. All of this is “general physics", not specific to 
phase field. It is the classical picture as drawn by Wulff [10], Herring [3], and others. 
But phase field incorporates this physics naturally. In the phase-field literature about 
dendritic growth, a simple anisotropy function with fourfold symmetry, in which the 
strength of the anisotropy is ô, is popular: 


c (8) = og (I + å cos [4(0 — 09))) (3.7) 
0*(0) = og I + 8 cos 4(0 — 09) — 165 cos 4(0 — 69)) 
= og (I — 158 cos 4(0 — 09). (3.8) 


Before ending this section, let us make several remarks: 

Firstly, we see that for ô > i the stiffness becomes negative (for 0 = 0). 
Phenomenologically speaking, inclinations of the interface that have a negative 
stiffness are forbidden. They collapse into edges and corners of the crystal: the 
crystal surface has “missing angles.” In phase fields, this is not possible to realize, 
since, due to the diffuseness of the interface, edges and corners are always rounded 
off with à maximum curvature of the order of the reciprocal interface width » 


Therefore, one has to "regularize" the corresponding function of the interface 
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stiffness with respect to orientation, cutting it to a positive minimum value of the 
stiffness. We regard this as "technical." 

Secondly, in 3D, i.e., having two inclination angles, one usually expands the 
interface energy and the stiffness in spherical harmonics. For a strongly anisotropic 
interface energy, e.g., for faceted crystals, other models may be appropriate (see 
Example at the end of this chapter). 

Thirdly, in a solid, we also need to consider the three misorientation angles for 
each inclination-angle landscape. 

Finally, at the atomistic scale, i.e., for nanocrystals, all these phenomenological 
relations—and also their microscopic or mesoscopic phase-field equivalents—will 
break down! The description of "capillarity" by curvature is only a statistical 
description of the atomic movement at interfaces. Atoms do not care about the 
curvature that we measure at an interface at a scale large compared to the atomic 
scale. But collectively they behave as if they would follow these rules. And there are 
other effects involved in the atomic motion connected to interfaces such as coupled- 
motion [1] or stress-strain effects at the nanoscale [2, 8]. Figure 3.1a schematically 
explains the meaning of the inclination angle with respect to the crystallographic 
axis and the normal vector of the interface in 2D. Figure 3.1b displays a simple 
model for the interface energy and stiffness as functions of the inclination angle. 


3.2.2 Interface-Energy Anisotropy: Phase-Field Picture 


In phase-field theory, the whole issue of interface-energy anisotropy is considered 
naturally, though there currently remain many open issues of technical character: 
it is not so easy. Let us start from the beginning, the phase-field equation as the 
functional derivative of the free energy, Eq. (2.6): 


So — 8F 


E x E (3.9) 


Here, we will only investigate the capillarity contribution, i.e., we will neglect the 
bulk energy difference between phases Ag = 0. In the physical notation, with 
anisotropic interface energy o (71) as a function of the normal vector to the interface 
n but with isotropic interface width 7, the variation of the free energy then becomes: 


2 
sr= f ax{ notin a e? + = ea - e] (3.10) 
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2 
+ nôch) [vor + - $d- 9| | (3.11) 


We use the double-obstacle potential; the case with the double-well potential 
works identically. In the first part (3.10), the variation is applied to the gradient 
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and potential contributions. This part we already know from Chap. 2, Eq. (2.15). 
The new part is the variational derivative of the interface energy (3.11) since it is 
anisotropic with respect to inclination, and the interface normal zi is a function of 


the gradient of the phase field n = S . First, we notice that the term in the square 


brackets of Eq. (3.11) is a positive function of $ and can be approximated, neglecting 
higher-order curvature terms and using relation (2.31): 


2 
[vor + m Ie — o»| ~ UVA). (3.12) 


A general discussion of the variation 50 (7) in 3D is quite involved, so we restrict 
ourselves here to 2D for the sake of tractability. Setting the MA angle to zero 
09 = 0, we define the angle of the surface normal 0 = arctan( 3> $y -) with the derivative 
of $ in the x and y directions of a Cartesian coordinate system. Then we expand: 


80 (0) (Vo)? = “2 save) (3.13) 
u do (0) ox dhy = $500 2 
= wer) (3.14) 
= 220 ($06, — $,6:) (3.15) 


O do (0) ðo (0) 
= n T nje- (5 $y) so) 


| 9?o(0) 98 
= ag (= bx D) bo 


_ 0a (6) 
~ 98? 


IVó| x 59. (3.16) 


We leave the proof of this relation as an exercise. To help with this: we have three 
different differentials to consider in expanding the variation ôo (0). Be reminded 
that all these differentials are linear operators that commute. In Eq. (3.13), the chain 
rule is used to expand the variation of o into a differential in 0 and a variation in 
0. Then, in Eq. (3.14), the variation in 0 is expanded in differentials in x and Øy, 
which by themselves are differentials in space. Then, partial integration is applied 
to the last term in Eq. (3.15) to respectively integrate the variations 6$, and d@y in 
space to separate the single variation 84 (again neglecting boundary integral terms). 
This needs a lot of bookkeeping, but it is more or less straightforward. The final step 
8 hy — a by = |V@|« is an approximation to lowest order in x, where the variation 
of the angle along one Cartesian direction is weighted by the phase-field gradient in 
the other direction. This variation, clearly, has to vanish for a planar interface where 
the curvature is 0. 
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Collecting all contributions, we end up with the desired expression for the 
evolution of the phase field (valid for 0 < $ < 1 and Ag 0), the equivalent 
of (2.36): 


2 
a = MP e von [vo > (å 2] “Ved Bach (3.17) 


ot 2 


Please note that throughout Chap. 2, the interface energy o was taken as isotropic, 
so no torque o” was considered. In general, this should always be there, as interfaces 
in crystalline phases are always dependent on orientation and inclination, even if 
only weakly. Several aspects shall be remarked upon. The form of Eq. (3.17) may 
not be the optimal choice for a numerical solution, since the prefactor o* = o + o” 
in front of the Laplacian can be strongly varying for strong anisotropies. If an 
explicit time-stepping scheme is used for the phase-field equation, this may lead to 
a significant reduction of the possible time steps and/or it may introduce instabilities 
in the front. An alternative is to evaluate the curvature of the interface locally and 
average it within a certain region. Then, one can use the Herring torque term o” as 
a smoothly varying driving force for the interface. Alternatively, one may directly 
discretize the variation óc (n), Eq. (3.11), as done, e.g., in [4]. A comparison of 
different approaches for dendritic growth is given in [5]. 


3.3 Exercises 


Exercise 

Repeat for yourself the relation of the mathematical curvature (3.1) with the 
corresponding phase-field expression (3.4). Consider that this relation is only 
true if the phase field is in the "right contour." 


Exercise 
Prove the derivation of the Herring torque from the phase-field functional 
(3.16). 
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Example—Equilibrium Shapes 

Anisotropy of interface energy is important in materials because it signifi- 
cantly impacts the morphology of the grain surface, the kinetics of interface 
migration, and the overall structure of the composite. According to the 
conventional definition of interface-energy anisotropy, it can be described 
as a function of the crystallographic orientation of the interface, which 
results in the equilibrium forms of individual crystals given by the Wulff 
construction [10]. We use an anisotropy model here, with the interface energy 
being a function of the inclination angle 0 only. The interface stiffness with 
respect to the inclination angle is modeled as [7]: 


aola? 


Ge = og (Og) + o" (Oy) = 


(3.18) 


uw 


NI 


(sin? (65) + a? cos? (04)) 


Validation of the Equilibrium Shapes 

The proposed model was employed in a 3D simulation box with 64? grid 
points of an initially spherical grain inserted into the melt. The inclination 
angle is computed using the interface normal of the phase field and different 
basic sets of facet normals, as shown in Fig. 3.2. 


© 9 @ @ 


a) Hexahedron b) Octahedron c) Cuboctahedron d) Rhombicuboctahedron 


Fig. 3.2 Equilibrium shapes with different facet vectors: (a) only {100} facets; (b) only {111} 
facets; (c) {100} and {111} facets; (d) {100}, {110}, and {111} facets 
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Further Reading 


* Rigorous derivation of the phase-field Herring torque: [6]. 
* Cahn-Hoffmann £-vector formalism: [9]. 
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Chapter 4 A) 
Temperature TEEM 


4.1 Thermodynamically Consistent Derivation of the 
Phase-Field Equation Coupled to Temperature 


In the previous chapters, we studied the phase-field equation itself; it can be 
used to propagate planar waves, but also to study the effect of capillarity on the 
equilibrium structures of crystals. For phase transformations, however, we have to 
consider the release/consumption of latent heat and solute redistribution in an alloy. 
In this chapter, we will study the interplay between the growth of a crystal and 
the temperature field around it. We will do this for a solidification problem, but 
the situation is not restricted to solidification. During a solid-state transformation, 
the release of latent heat also has a significant effect. We will do this for a pure 
substance in which solute redistribution plays no role. However, in general, both 
fields—temperature and solute (see Chap. 5)—have to be coupled to the evolution 
of the phases in an alloy. See also “Example—Temperature Evolution During Rapid 
Solidification" at the end of this chapter. 

The Gibbs free energy in the case that we solve the heat-conduction equation 
is not an appropriate starting point, since the Gibbs energy requires that the 
temperature be given. Instead, one starts from the entropy functional S. Following 
Wang et al. [5], the entropy functional S as the integral of the entropy density s 
over the domain €2 is defined by the internal energy density e and the free-energy 


density f: 
s= fs= | (4.1) 
2 o T 
e — pcpT + L(1— 9). (4.2) 
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42 4 Temperature 


The internal energy of a solid is assumed to be linearly dependent on temperature 
T with constant density p, specific heat capacity cp, and latent heat of fusion L. 
The governing equations for T and $ are derived consistently with the principle of 
entropy production $20 


" 5 f 1 [af af 
= a (he) TE Vane. s 
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Inserting the energy model (4.2) into (4.4), we obtain: 


: Mr 
é = pep T — Lå = VZ VT 
pcpT = VArVT + Lø, (4.5) 


which is the well-known heat-conduction equation, with thermal conductivity Ar — 
A The phase-field Eq. (4.3) is then almost identical to the previously derived 
equation, e.g., for the double-obstacle potential (2.36). The only modification is the 
prefactor +$, which can be assimilated into the phase-field mobility M?. The ther- 
modynamically consistent derivation of phase-field models is of special importance, 
because it enables the correlation of the model parameters with each other, as well 
as the establishment of a sound theoretical background in thermodynamics. 

The derivation above defines a pair of coupled partial differential equations. Both 
of these are parabolic diffusion-type equations. The interesting feature is the source 
term in each equation: the release of latent heat in the heat-diffusion equation, which 
depends on $, and the driving term for the phase field, which depends on the actual 
temperature of the interface T. This mutual coupling leads to the morphologically 
unstable evolution of a dendritic solid-liquid interface. It also poses a challenge 
for numerical solution schemes of this set of equations. Usually, they are solved 
explicitly in a staggered scheme, i.e., sequentially: one and then the other. In the 
literature, more advanced schemes can be found, but these do not prevail in general 
practice. A simple piece of code in C++, a typical 100-line program (not counting 
input and output) for a phase-field problem, is given in Appendix A.1 for you to try 
yourself. 


4.2 Thin-Interface Limit 


Let us now focus on one specific problem associated with a physically meaningful 
solution of the above-outlined solidification problem: solidification occurring at 
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relatively high temperatures with high mobility of the solid-liquid interface Må. 
Typically, this transformation is treated as "diffusion controlled,” meaning that 
the interface T; is set to capillarity-corrected thermal equilibrium at temperature 
Tm + ist in which x < 0 for a convex dendrite tip and A Ss, is the entropy of 
fusion. This condition tells us that the Gibbs- Thomson equation, which we have 
used to control the speed of the interface, becomes undetermined. As a reminder, 
the Gibbs-Thomson equation reads: 


v = Må, (o*« + ASs (Ti — Tm)) - (4.6) 


For any finite velocity v, the condition Mg, — oo determines only the interface 
temperature Ti = Tm — Be. but the velocity cannot be specified. To proceed in 
this situation, we consider the temperature profile ahead of a growing solidification 
front, as sketched in Fig. 4.1. 

In the growing solid on the left-hand side of Fig. 4.1, the temperature can be 
assumed to be uniform. In the sharp-interface picture (blue line), there is a kink in 
the temperature profile, which relates to the release of heat at the solidification front. 
In this case, it helps that the kink of the temperature at the front can be related to the 
interface velocity (see the classical Stefan problem [4]): 


oT oT oT 


Lv = As | | 
Ox IS Ox IL Ox 


- (4.7) 
in which |s and |r, denote that the temperature gradient is evaluated in the solid 
or liquid, respectively. Fluxes in the solid can be neglected under steady-state 
conditions because the solid is uniformly at the melting temperature. Equation (4.7) 
is a balance equation between heat release due to solidification Lv and heat 
extraction due to diffusion into the supercooled liquid Mår. We can see that the 
velocity is now calculated from heat diffusion instead of being proportional to a 
thermodynamic driving force (Gibbs-Thomson condition with finite mobility). This 
is called the “diffusion-controlled limit.” 

In phase field, this limit seems difficult, because the phase-field equation for $ 
corresponds to the Gibbs-Thomson equation. The seminal contribution of Karma 
and Rappel [1, 2] was to realize that the temperature-diffusion equation, or heat- 
conduction equation, can “easily” be integrated into the phase-field equation! 

To understand this, we go back to Fig. 4.1: the dotted line indicates the phase- 
field contour, which is diffused over the interface width n; the green line indicates 
the temperature corresponding to this solution. The release of latent heat is no 
longer concentrated at the front (producing the kink in the temperature profile in the 
sharp-interface picture) but is continuous over the interface. Correspondingly, the 
temperature profile is smooth. However, outside the interface, we demand that both 
temperatures, the sharp-interface and diffuse-interface temperatures, match. This is 
called the “asymptotic matching condition,’ as applied to solving the temperature 
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Fig. 4.1 Sketch of the temperature profile through an interface of a pure substance growing into 
an undercooled melt. The sharp front position is located at the center of the interface with a kink 
in the sharp-interface temperature. To the left is the solid at melting temperature Tm (neglecting 
capillarity effects for the moment); to the right is the undercooled melt, and the front is moving 
from left to right. The solid blue line depicts the temperature of the sharp-interface problem with 
a kink due to the release of latent heat at the moving front. The dotted line indicates the phase 
field, which is smeared out over the width 7. The green line corresponds to the temperature profile 
of the phase-field model with a smeared-out release of latent heat within the interface. This shall 
match the sharp-interface solution outside the interface, but it has a systematic deviation within the 
interface. The latter produces a “spurious undercooling;" which shall be used to compensate for 
"numerical undercooling" (see main text) 


profile for a given phase-field solution analytically (in a 1D direction normal to the 
interface). 

We can see directly that, systematically, the temperature of the phase-field model 
has to be lower (in the case of solidification into an undercooled melt) than the 
respective sharp-interface temperature. On the one hand, this is not good because it 
is a systematic error that cannot be avoided, and errors are never good; on the other 
hand, as long as the temperature outside the interface is correct, we will accept it. 
This is consistent with the notion that a mesoscopic phase-field model does not 
claim full physical resolution inside the interface. This systematic error then turns 
out to be very good from a practical point of view: it helps us to relate the velocity 
of the interface to undercooling! 
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There is no physical undercooling of the interface in the diffusion-controlled 
limit of a sharp-interface model, but there is necessarily a spurious undercooling 
in the case of a diffuse-interface model. If the right temperature solution as 
function of the velocity can be calculated analytically, it should coincide with the 
spurious undercooling of a good numerical solution Ag""mericál For the analytical 
solution, please refer to the literature [1—3]. Here, we simply state that the spurious 
undercooling can be (on average over the interface) expanded to lowest order in 
the front velocity v as Ag" ~ Av with a constant A depending only on the 
materials parameter Å, ASsr, and 7. This can be understood easily because the 
spurious undercooling must vanish if the system reaches equilibrium, i.e., v — 0. 
With little math, we can reformulate the Gibbs-Thomson equation: 


v = M? [0*k + Ag] (4.8) 


M? [o*k + Ag + 0] 
M 


$ ok +Ag+ Ap red EE Agteurieus] . 


Here, again, we add an “intelligent 0", 0 = Ag""merical — A gspurious We assume 
that our numerics are good and that the numerical solution of the “thin interface 
temperature" within the interface matches the correct analytical solution (green 
line in Fig. 4.1). Then the numerical driving force can be read from the numerical 
solution and the spurious driving force can be handled analytically as a "spurious 
velocity". This we do! 


v(1-F M? A) = M? oe + Ag + V i 


M? * numerical 


This defines the effective mobility of a diffuse-interface model, M ia = De 
It is largely determined from diffusion in the dying phase. This means that a phase- 
field model, coupled to temperature diffusion, runs with a spurious undercooling and 
a finite effective mobility, even if the physical mobility M? — oo. It also matches 
the sharp-interface solution outside the “thin” interface. This is called “thin-interface 
asymptotic.” The driving force Ag* comprises the physical and the spurious (but 


necessary) numerical contributions. 
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4.3 Exercises 


Exercise 
Repeat (4.8) to (4.9) for yourself! 


Example—Temperature Evolution During Rapid Solidification 

In this example, the system is subjected to the solidification conditions of an 
additive manufacturing process, i.e., a process with a very high cooling rate 
and very high thermal gradients. In the coupled phase-field model, the heat 
conduction equation, 


pCyT = V - (AVT) + Qrocai, (4.10) 


is solved implicitly. Here, Qoca] handles all kinds of heat sources, including 
both the release of latent heat and external heat sources. From the simulation 
results shown in Fig. 4.2, we can see the growth of dendrites in cellular form 
owing to the very high temperature gradient and cooling rate. The line scan 
of the temperature plot along the z axis highlights that the dendrite tip is the 
hottest region during solidification. This is due to the release of latent heat 
during phase transformation. 
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Fig. 4.2 Phase-field simulation results obtained for a solidification process under additive man- 
ufacturing conditions. Left: evolution of solid phase. Middle: temperature distribution at the box 
surface in color coding for a given time step. Right: line scan of temperature along the z axis with 
a peak temperature at the dendrite tip 
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Heat Release During Solidification 


Further Reading 


* The original publication on the thin-interface limit is by Alain Karma and 
Wouter-Jan Rappel [2]. The limit for a phase field coupled to solutal dif- 
fusion (see Chap. 5) is presented in [1]. Here, the so-called "anti-trapping 
current" for a model with a strong difference in the diffusion coefficients is 
also introduced. In alloy solidification, the liquid diffusivity is much larger 
than the solid diffusivity, and this has important consequences for treating 
the thin-interface limit. 


References 


1. A. Karma, Phase-field formulation for quantitative modeling of alloy solidification. Phys. Rev. 
Lett. 87(11), 115701 (2001) 

2. A. Karma, W.-J. Rappel, Quantitative phase-field modelling of dendritic growth in two and three 
dimensions. Phys. Rev. E 57, 4323-4349 (1998) 

3. I. Steinbach, Phase-field models in materials science. Model. Simul. Mater. Sci. Eng. 17, 073001 
(2009) 

4. C. Vuik, Some historical notes on the Stefan problem. Nieuw Archief voor Wiskunde IV 11(2), 
157—167 (1993). ISSN: 0028-9825 

5. S.-L.Wang et al., Thermodynamically-consistent phase-field models for solidification. Physica 
D 69, 189—200 (1993) 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate 
credit to the original author(s) and the source, provide a link to the Creative Commons license and 
indicate if changes were made. 

The images or other third party material in this chapter are included in the chapter's Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not 
included in the chapter's Creative Commons license and your intended use is not permitted by 
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. 


Chapter 5 A 
Concentration FECA 


The last chapter discussed the coupling of the phase field to temperature. You will 
remember that the phase-field functional has three contributions: the “gradient,” 
the “potential,” and the “bulk free-energy difference.” The gradient has to fight 
with the other two contributions: the fight with the potential represents capillarity, 
and the fight of capillarity with the bulk free-energy difference determines phase 
transformations and growth or shrinkage of phases. Therefore, all these ingredients 
have to be considered simultaneously. 

In this chapter, we consider the effect of solute on phase transformations. The 
treatment of solute within the diffuse interface of a phase-field model is one of 
the most controversial issues in phase-field theory. Here, we will give a condensed 
statement of the problem and then elaborate on the most recent model to couple 
phase-field theory and solute diffusion within the mesoscopic picture, the so-called 
“finite interface dissipation” (FID) model [12, 14]. 


5.1 Phase Field and Solute Concentration 


The problem with treating the phase field and solute concentration in a consis- 
tent theory relates to the physical question of how crystal structure and solute 
composition correlate during a phase transformation. This is the chicken and egg 
problem: who comes first? During a phase transformation in an alloy, the crystal 
structure typically changes along with the solubility of the alloying elements. 
We can likely assume that both mechanisms go hand in hand, but a generally 
accepted rule cannot be given here. A special case is “spinodal decomposition,” 
in which the crystal structure stays the same, but due to a miscibility gap, the 
solute partitions into two composition sets. The theory for the chemical interface 
energy that dominates in this case was proposed by Cahn and Hilliard [1]. This is 
now frequently applied in microscopic phase-field models of phase transformations 
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with structural components. This may work well if: (1) the structural part of the 
transformation is considered together with the chemical interface energy; and (ii) the 
interface is treated at the atomistic scale. 

In a mesoscopic phase-field model, in which the interface width may reach the 
micrometer scale, a Cahn-Hilliard approach to solute diffusion will cause severe 
errors because the chemical interface energy in a Cahn-Hilliard model scales 
with the interface width [11]. If we consider the interface at the atomistic scale, 
with a roughly 1-nm width compared to a scaled interface width of 1 jum in a 
mesoscopic simulation, this error will be a factor of 1000. This situation is depicted 
in Fig. 5.1: the Gibbs energy of a two-phase system f®°™ has a hump in the two- 
phase region between the equilibrium compositions c£ and c$ if the composition is 
treated as a continuous variable c between c$, and ch- The integral over this hump 
in the normal direction through the interface gives the chemical interface energy 
Ochem = f dn je xn. 

There are three options to cope with this problem: 


e Use adaptive methods to scale down the interface to the atomistic scale. 

e Construct a Gibbs-energy landscape without a hump between the phase concen- 
trations. 

* Split the composition c into phase compositions cy. 


Fig. 5.1 Sketch of the Gibbs energy of a two-phase system (left). If we map the concentration 
between the phase concentrations cg and cg onto the normal n through the interface, the hump in 
the chemical Gibbs energy f *'* transfers to a region in space with increased energy (right). This 
region scales with the interface width 7 
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The first approach, to use adaptive finite elements with mesh refinement towards 
the interface, raised some interest in the early days of the application of phase 
fields to solidification problems (see, e.g., [10]). This works sufficiently well in 2D; 
however, it is hardly possible to resolve structures in the micrometer range in 3D. 

The second approach, as applied by Echebarria et al. to dendritic solidification 
of a binary model alloy [2], is difficult to generalize for an arbitrary Gibbs-energy 
landscape, and particularly for multicomponent materials. Nevertheless, it presents 
a rigorous procedure to get rid of spurious effects of coupling between solute 
redistribution and phase evolution in the interface. 

The third approach is mostly used today for applications in multicomponent 
materials. First introduced by Tiaden et al. [13], this splits the overall composition 
c into the phase compositions cy and cg, where $ indexes the o phase: 


c= bly + (1 — $)cg. (5.1) 


Hereby, the phase compositions cy and cg may deviate from the compositions in 
equilibrium c and ch- The extra degree of freedom that arises from this construction 
is fixed by relating the phase compositions by a partition coefficient taken from 
the equilibrium phase diagram, as done in the original work [13]. A generalization 
with an extrapolation scheme to multicomponent phase diagrams is presented 
by Eiken et al. [3]. An alternative is to postulate local chemical equilibrium 
within the interface, as proposed by Kim et al. [5], i.e., equal chemical potential: 
Ma = Up = p. This is called the “KKS” model from the names of the authors, and 
it is quite popular. It has specific problems of a technical nature when applied to a 
general Gibbs-energy landscape. Therefore, Plapp developed the so-called grand- 
potential approach [8]. 

The grand potential is a Legendre transformation of the free energy, in which 
the chemical potential replaces the composition as a state variable. Now, if we 
assume equal chemical potential within the interface, there is only one variable, 
u, and we are automatically in the minimum state of the free energy within the 
interface. Thereby, as a side product, the spurious interface energy is gone. Now 
one solves a diffusion equation for the chemical potential (instead of the diffusion 
equation for the composition) together with the phase-field equation. The diffusion 
equation for the chemical potential is formally identical to the diffusion equation 
for the temperature [Chap. 4, Eq. (4.5)] with a source term proportional to the rate 
change of the phase field $. One difficulty remains: how to invert the chemical 
potential into the phase compositions, which are, of course, the variables of interest 
if we are investigating a phase transformation in an alloy. An oft-applied remedy 
is to approximate the Gibbs energies of the individual phases as parabolic. The 
chemical potential then becomes linear in concentration, and one can directly read 
the composition in the bulk phases from the chemical potentials. 

One severe drawback of all the approaches described above is that they are 
restricted to local equilibrium within the interface, or at least to “close to local equi- 
librium.” If this condition breaks down—e.g., in rapid solidification, where solute 
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trapping appears, or in solid-state phase transformations at lower temperatures (such 
as bainitic or martensitic transformations)—one has to take a different approach. 


5.2 Finite Interface Dissipation Model 


The general idea of the FID model [12, 14] is to treat the phase compositions 
as separate variables within the interface that are connected by the conservation 
constraint of composition. Technically, one separates the diffusion fluxes into 
long-range fluxes within the bulk phases and into short-range redistribution fluxes 
between the individual phases within the interface. Then one can start from arbitrary 
initial phase-composition conditions; they shall converge toward local equilibrium 
within the interface if the kinetics of the process permits. We also assume that, in 
reality, the composition within the interface is not strictly in local equilibrium. 
To build such a model, we start from the following assumptions: 


e The phase compositions c, Within the interface are a constraint to the mixture 
composition c = 5^, daca. For a solid-liquid interface, this is simply c = 
Pscs + ØLCL 

* The phase concentrations are treated as independent variables, only subject to the 
sum constraint above. 


In a mesoscopic model, the interface is taken as an effective reference volume 
without specifying its actual position and orientation. Within this reference volume, 
we take the two phases as “mixed” and allow them to exchange solute to lower 
the total Gibbs energy of the reference volume. We define the chemical free-energy 
density /*"*" as a weighted sum of the chemical free energies of the individual 

chem. 


phases f7 Ip 


FO” = ha fa ™ + dp fg" +A (c — Paca + beep) ined 


where the Lagrange multiplicator Å ensures conservation of the concentration and 
gives us the possibility to take the phase concentrations cy fully as independent 
variables. Then, we define a redistribution flux between the phases from the demand 
of energy reduction with the permeability P (for details see [12]): 


8 à à 
duce epo pq p ba — bya), (5.3) 
Co Ocq 3Co 
. ô 3 fp 
-— ff Sap hem SP go he | 54 
Ppcp cg cy f |o T op | (5.4) 


The permeability is the inverse of the resistivity of the interface against redistri- 
bution between the different phases. We can certainly assume that the interface has 
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an influence on the diffusion of solutes, that it hinders diffusion or even makes it 
easier because of free volume within the interface region. This will depend on the 
particular kind of interface. For a high permeability, it is easy to exchange solute 
between the phases. A low permeability will characterize a passivated interface 
through which no solute can pass. 

From the conservation condition ¢ = 0, which shall hold within one reference 
volume without external fluxes (we will add these later), one determines the 
Lagrange multiplier: 


0 = (paca + bpcp) = baca + bpcp + Pala + Pptp 
= bala + dgcp + bala 
0 
18 ined dpa] l (5.5) 
CB 


ð 
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ð 
i= paz (5.6) 
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Substituting A into Eqs. (5.3) and (5.4), we get the evolution equations for each 
phase inside the reference volume, 


; 0 a N 
dala = Papp E: = am) T Pa Pa (cp — Ca), (5.7) 
0cg Cy 
: a Of, : 
deg = Poupa ( La — 8) + øpdglen = cg). 5.8) 
Ocq deg 
Finally, we add the long-range diffusion in the bulk phases in the standard way. Da 
and Dg are the diffusion coefficients, and Lg = ye and ug = ab are the chemical 
potentials: 
Pata = V (ha Da Vea) F Papp (us — Ba) + $a. (cg — Ca); (5.9) 
øpåp = V(Óg Dg Veg) + Padg(V — Vg) + $pbp (ca — cg). (5.10) 


By construction, both equations overlap in the interface, where two mechanisms 

of solute redistribution are considered: redistribution due to a chemical-potential 
difference pe — pg Å 0 and renormalization if the fraction of phases changes 
$a = —óg # 0. 
Remark All phase-field implementations that use the splitting of the composition 
c into phase compositions cy have to include a redistribution step and a renor- 
malization step. This redistribution procedure is, unfortunately, never discussed in 
detail in the literature, but it is considered as a technical concern. Partitioning and 
redistribution is, however, a crucial consideration if you want to implement your 
own solution. 


54 5 Concentration 


The diffusion Eqs. (5.9) and (5.10) for the phase compositions œ and f converge 
to the KKS model in the case of a very high permeability P — oo because the 
redistribution is then very fast and the chemical potentials jj, and ug must become 
equal. In the case of a low permeability, however, the interface concentration and 
the partitioning in the interface may significantly deviate from local equilibrium. 
Currently, there is no general model to determine the permeability parameter for 
a particular interface. The model as presented in the original paper [12] must be 
considered to be in error, since it has an inverse dependence of the permeability 
on the interface width, which must be considered as nonphysical in a mesoscopic 
model. A special model for the permeability in an electrochemical system, derived 
from experiments, is given in [9]. 

The model has another intriguing feature regarding the influence of the perme- 
ability on the motion of the interface. Going back to the model for the chemical 
free energy (5.2), we notice first that the Lagrange parameter A takes the role of 
an effective chemical potential of the interface. In fact, it is—according to (5.6)—a 
mixture of the chemical potentials of the individual phases, linear in $. Furthermore, 
it is also a function of the rate change of the phase field $! Since all kinetic-model 
equations—phase-field, temperature, and composition up to now (stress and strain 
will be added in Chap. 2)—are derived as thermodynamically consistent from the 
same free-energy functional F, this term will enter all kinetic equations. The phase- 
field equation [cf. Chap. 2, Eq. (2.36)] becomes: 


; ôF 
= —mM? 
$a = Mg 
? 1 
=M9 [ov - 7 G -ØI-TV$U-Øag |. GID 
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Agg = fa — fg — Ma — cp) 


Eg (ov. +(L-ø)ug Ma <2) (ca — cg). (5.12) 


Rearrangement of the terms proportional to ġa finally leads to: 
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The phase-field mobility is modified by a mechanism that is related to solute 
redistribution and diffusion. The interface mobility M? in the thin-interface limit 
should be taken as the effective mobility derived in Chap. 4, at least for high 
permeability P when the model converges to the KKS model. To date, a consistent 
treatment for the non-equilibrium case (small P) for thin interfaces is still missing. 

We see that there are three distinct limits for the new phase-field mobility K: 


. P — oo: K —^ M9; 
* P —> 0U (cg — ca) #0: K > 0; 
* P —> 0U (cg — ca) = 0: K > M?. 


The first case, with high permeability, leaves the phase-field equation unchanged, 
and the ID model converges to the KKS model with equal chemical potentials within 
the interface. The second case represents a passivated interface: no redistribution of 
solute is allowed, but there is a jump in concentration between the phases. Any 
transformation is forbidden, and the transformation stops. The most interesting case 
is the third case: there is no concentration difference between both phases, as in 
a partitionless transformation like martensitic transformation, and the phase-field 
mobility is unaffected, i.e., the transformation can proceed with very high speed. 
Again, no detailed studies of this have been conducted to date. 


5.3 Multicomponent Alloy Transformation 


Up to here we were mainly speaking about a binary alloy with one composition 
variable c. This is important and gives you all the relevant information about “alloy 
transformation,” including the general problems associated with it as partitioning 
and getting the phase diagram right. The extension to multicomponent alloys is 
straightforward (see [14]) from a theoretical point of view. It is more than “involved” 
from a practical point of view. Technical alloys, steels, nickel based superalloys, 
aluminum or copper alloys, brass and others, but also minerals and ceramics are 
multicomponent with 10 or more individual components. Also vacancies have to be 
considered for completeness. The first challenge which arises is the phase diagram 
from which we would like to read the equilibrium compositions of two particular 
phases, see Fig. 1.1: There is no way to display a phase diagram of a 10 component 
alloy and we need a numerical procedure to do the calculation of equilibria and 
deviation from equilibrium automatically. 

The second challenge is related to redistribution and diffusion during a phase 
transformation. The individual components i of an n,-component alloy in phase 
state o, are not independent but connected by constraint to the sum constraint: 


ne 


waa (5.16) 
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The third challenge, and a challenge in its own right, is the construction of the 
Gibbs energy landscape of a multicomponent alloy with 10s of possible phases. 
For the latter we have to rely on the so-called CALPHAD (CALculation of PHAse 
Diagrams) technique (see e.g. Lukas et al. [7]). Open or commercial databases are 
available which provide the Gibbs energy of a special phase (stable, metastable or 
even unstable) as function of temperature, pressure and all available components 
in a special machine readable format. These you simply insert into the phase-field 
functional (1.4) to calculate the Gibbs energy difference Ag($, T, c, €, ...). In the 
same way one can access the chemical potentials pg, Mg in the Sect. 5.2. Now 
they are also multi-component tg > TU and one has to take the derivative in 
the directions i [14]. As said, in practice this may become very involved. 

The last challenge to be mentioned here is multicomponent diffusion, again a 
challenge in its own right. We would like to incorporate atomic mobilities from 
so-called “mobility databases” which accompany thermodynamic databases. We 
reference here to "Further reading." 

In the literature several approaches to cope with “multicomponent transforma- 
tions" are established: 


* Hard-coding of thermodynamic functions. This will be helpful if you have a low 
ranked problem, say n < 4, and a suitable analytic description of the Gibbs 
energy landscape available. Also you will need your own code and sufficient 
programming skills. And you will have to do the job for each new problem. 

* Local parabolic approximation of the Gibbs energy landscape. Local means here 
in composition space for your given alloy system, temperature and pressure. The 
approximation you may do with the help of a CALPHAD software, or from 
a given analytic description. The approach has also the advantage that for this 
quadratic approximation the chemical potential is linear with the composition! 

* Local linearization of the phase diagram [3]. This approach is somehow similar 
to the previous approach, but quite different in technical aspects. We consider it 
the most efficient approach for large scale computation, but it may not be suited 
for larger deviation from local equilibrium at the interface. 

* Direct coupling to databases. This is the most general way of incorporating the 
full CALPHAD scheme into a phase-field calculation. Of course it has the highest 
computational cost. An example using the OpenPhase software is shown at the 
end of this chapter (Fig. 5.2). 
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Fig. 5.2 Phase-field simulation of directional solidification of a four component Ni-Al-Cr-Ta 
superalloy under additive manufacturing conditions. Left: Al composition on the surface of the 
simulation box. Right: Composition profiles for the individual alloy component on a cross section 
horizontally through the simulation box (position indicated in left figure). Colour bars above each 
profile display the composition range. The directed dendrites under these conditions show almost 
no side branches and interdendritic y’ is precipitated at the end of solidification 


5.4 Exercises 


Exercise 
Derive the expression for the Lagrange multiplier A (5.6). 


Example—Multicomponent Alloy Solidification 

The following example relates to multicomponent alloy solidification in a Ni- 
based superalloy. A model alloy of the commercial alloy CMSX4 with four 
components, Ni, Al, Cr and Ta, is constructed such to yield a comparable 
fraction of secondary y’ phases at the solvus temperature of this alloy. The 
solidification conditions relate to additive manufacturing, as in the previous 
Example in Chap. 4. 
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5 Concentration 


Multicomponent Alloy Solidification 


Further Reading 


* Extrapolation scheme for multicomponent phase equilibria [3]. 

* Pair-wise multicomponent diffusion approach [4, 6] 

* FID model for multicomponent alloy transformation [14] and sublattice 
ordering [15]. 
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Chapter 6 A 
Multi-Phase-Field Approach erts 


6.1 Phase-Field Functional for Multiple Phases 


Technical materials have multiple phases. In a CALPHAD (CALculation of PHAse 
Diagrams) database for iron alloys and steel, there are Gibbs-energy functions of 
several hundreds of different physical phases tabulated, and these are dependent on 
the crystal structure, composition, temperature, and pressure. According to Gibbs’ 
phase rule, there can, however, only be n = 2 + n, stable phases in equilibrium, 
where n, is the number of components of the alloy composition. For a steel, we 
may put n, = 10, i.e., 12 different phases could be stable in equilibrium. Since, 
however, a steel is rarely in thermodynamic equilibrium but is more commonly in a 
kinetically stabilized off-equilibrium state with a locally varying composition, you 
may find 20 or more individual phases in a steel sample. In phase-field theory, we 
also attribute different grains—i.e., crystallites with the same crystallographic phase 
but a different orientation in space—to different phase fields. We do this to identify 
the interfaces between these grains, which will evolve in time (see Chap. 3). We can 
then have thousands of phase fields. How do we tackle this? 

For completeness, there are two fundamentally different approaches for this 
"multi-grain" problem. The microscopic, or physical order-parameter models, use 
a free-energy landscape. Generally speaking, this is a special potential function 
with different minima that are dependent, for example, on the orientation of the 
crystallites [1, 11]. 

Another approach is the introduction of orientation as an order parameter [5, 6]. 
These approaches are generally considered as elegant from a theoretical point of 
view, but they are difficult to apply to practical problems. The reader should build 
his/her own opinion on this. 

The third approach, which is more “brute force,” is to address each crystallite, 
regardless of its phase, by its own phase field. We call this the “multi-phase-field 
approach" [12, 13]. Also popular is the so-called *vector order parameter" model 
by Fan and Chen [2]. Both models consider a set of N phase fields ġo (a = 1... N) 
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Fig. 6.1 Scheme of a triple junction between three phase fields. Within the bulk regions, only 
the information of the phase-field index has to be stored to identify the phase. Within the diffuse 
interface, two phase indices and one value of the phase field have to be stored. In junctions, we 
need the information of all phase-field indices, and all phase fields need to be stored. Of course, all 
of this information evolves dynamically during the calculation when the interfaces move 


that can be attributed to different properties such as orientation, composition, crystal 
structure etc. We will describe here only the multi-phase-field approach [12, 13]. We 
distinguish three cases, as depicted in Fig. 6.1: 


* Bulk: only one phase field has the value of 1, all other fields have the value of 0. 

* Dual interface: two phase fields have values between 0 and 1, all other fields 
are 0. 

* Junctions: we define a number N < N of non-zero fields that overlap in a 
multiple junction. In Fig. 6.1 Ñ = 3 in the triple junction in the center. 


Within this lecture, we will only deal with "junctions," since the rest is known 
from dual phase fields as treated in the previous lectures. There is nonetheless 
enough to do... 

The first thing to do is to define the free-energy density of a multi-phase-field 
model: 


N N 
8 og å 
=D) X a [wvevesztee]e dase}. 6p 
a=1 | B=a+1 


We sum over all N phases/grains present within the junction for both the capillarity 
term, or interface, and the bulk free energy FU, weighted by the local fraction 
of this phase o given by the phase field øw. This bulk free-energy density will 
be a function of temperature, concentration, stress, and strain, and it may contain 
magnetic or other contributions [cf. also Eq. (2.9) in Chap. 2]. 

The capillarity term is more complicated; it is expanded in pairs of phases, i.e., 
we have a second sum over all N phases. Note that the second sum runs only from 
B — a 4- 1... N so as not to count contributions twice. There is also no diagonal 


term 6 = o because this form does not relate to an interface between two phases. 
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For triple junctions between three phase fields (usually a "triple line" in 3D), we 
find three pairs of phases. For quadruple junctions between four phases, we already 
have six pairs, and there are ten pairs for N = 5 and so on. The junctions are 
then modeled by the superposition of dual interface contributions with the interface 
energy of a pair ogg. In general, we have N combinatorial 2 interfaces for V phases. 

There may be higher-order contributions, i.e., contributions where more than two 
phase fields and their gradients are considered simultaneously. They should give 
a special penalty to multiple junctions, as discussed by Miyoshi and Takaki [7]. 
The present form (6.1) is just a pragmatic setting whereby all parameters are, in 
principle, well-defined physical entities. The interface width n is not considered a 
physical entity at the mesoscopic scale and is set to a constant value. The interface 
density [7-2 Vo. Vóg + 1 bad I] is expanded in pairs as a pragmatic choice. One 
can easily check (and please do so) that this form reduces in the case of a dual 
interface Ñ = 2 to the standard form with a double-obstacle potential (2.35) where 
the bulk free energy is weighted linearly in ġa. 


6.2 Double Obstacle versus Double Well in Multi Phase Field 


You may skip this part on first reading, but it is important! Why do we choose a 
double-obstacle potential instead of a double-well potential? An important feature of 


the multi-phase-field theory [12] is the sum constraint in junctions: E 1 9o = 1. 
For the double-well potential we write the potential (without prefactors) in 


o= P» Hab (6.2) 


a=1...N,B=a+1...N 


The maximum of the potential in the center of the junction, where 


1 
$^ = of — * for all a, B, ..., is then evaluated: 
N\ 1 1 = 
DW = WE x G for N > 1. (6.3) 


This means that for N > 3, the energy of the junction decreases with the order 
N and approaches 0 for large N. This must be termed “unphysical.” since junctions 
between objects lose their penalty, and the system would return to the disordered 
state. The double-obstacle potential is introduced [12] to remedy this problem: 


feo = D adl. (6.4) 


a=1...N,B=a+1...N 
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This has the same topology as the double-well potential (see Fig.2.2) but a 
maximum power of two. We calculate the maximum potential of the junction Tp 


29 = (S) «1 for N > 1. (6.5) 


This means that the energy of the junction increases with the order N and 
approaches a constant for large N, as it should. The main drawback of this potential 
is the non-analytical form with the absolute signs. But this is technical... 


6.3 Multi-Phase-Field Equation 


The most important feature of the multi-phase-field theory is the consideration of a 
sum constraint: 


Y ha =l. (6.6) 


The system must be closed in itself so that “holes” do not arise, and this has 
important consequences for further model development. The phase fields within a 
junction cannot be varied independently; m 0. Starting with a relaxation ansatz 
(2.6), we must apply the chain rule: 


Ibu ô dpp à 
ap ita * 2^ 06,36; F. (6.7) 


IDE: ' its 
The variation å. however, is generally unknown. Therefore, so-called "interface 


fields” Yag are introduced [12]: 


ô ô 


yee | eee E 6.8 
Wap P sl (6.8) 


These have the important property of automatically conserving the sum constraint 
(6.6) (for details see [3, 12]). With the phase-field mobility M$; defined for a pair 
of phases, we reformulate (6.7): 


$ ó 
Iba Mop i ô 5 | 
=-V E p=- F. 6.9 
ar 4g OUT ET e 


Now we have consequently decoupled a multi-phase problem into pairwise 
contributions. We end this lecture with two remarks, or warnings: 
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aL Ja 


(a) (b) 


Fig. 6.2 Force balance of a quadruple junction between four grains with isotropic interface energy. 
(a) The benchmark test starts with a simple configuration of rectangular grains. The triple lines 
between three phases are displayed. The inset shows the phase configuration in color coding. (b) By 
reduction of interface energy, the system relaxes to the final configuration with equal angles of 109? 
between the triple lines and a symmetric grain structure 


* Aselegant as it looks, it is tedious to implement in phase-field software: since F 
is itself expanded in pairs, there is a third loop over phases y 4 a  B. 

* As tedious it is to implement, it is nonetheless crucial if you are considering a 
real materials problem. Taking any shortcuts to avoid the general implementation 
will lead to a deadlock regarding the solution of the materials problem. 


A kind of benchmark for the equilibrium configuration of junctions is displayed 
in Fig. 6.2 where a brick-like initial structure relaxes into the equilibrium configura- 
tion [3, 4]. Note: This equilibrium configuration does not fit nicely into a rectangular 
box. Therefore one has to employ some tricks regarding boundary conditions and 
keeping the center within the box, it has the tendency to drift out, of course to 
reduce interface energy. Such important details for practical applications, however, 
are seldom discussed in close detail in the literature. 


6.4 Exercise 


Exercise 

Prove that the form (6.8) for the interface field automatically conserves the 
sum constraint (6.6) by adding a Lagrange multiplier A ibat Pa — 1} to the 
free energy F, see [12]. 
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Example: Anisotropic Grain Growth 

A multi-phase-field 3D grain growth model is formulated in [10] to simulate 
grain growth with anisotropic interface energy. The interface energy is defined 
as a function of the inclination angle of the boundary plane. It is found that 
the resulting grain-growth kinetics are strongly influenced by the anisotropy 
of the interface energy. This results in a slower growth rate and a distinct 
morphological evolution, as indicated by the presence of more cubic grains 
and a more significant number of triple-junction angles between 90? and 
180°, as illustrated in Fig.6.3. Additionally, the model closely matches 
various experimental results for NaCl and MgO polycrystalline minerals, 
where the distribution of grain-boundary planes peaks at low-index (100]- 
type boundaries. 
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Fig. 6.3 Three dimensional simulation of anisotropic grain growth in 3D [10]. A 2D cross-section 
and the bow up of an insert shows clearly the deviation of the angles of interfaces at junctions from 


the isotropic angle expected for isotropic grain growth. Also the structure shows a pronounced 
cubic texture 
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Further Reading 


* Amulti-phase-field model without the sum constraint, but with a consistent 
treatment of the bulk-energy part, is presented in [8, 9]. 
* Force balance at junctions [3, 14]. 
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Chapter 7 A 
Stress-Strain and Fluid Flow FEEN 


7.1 Coupling to Stress and Strain 


In this last lecture, before we change gears to the “quantum-phase-field” model, 
we want to (almost) complete the picture of the effects to be considered in real 
materials. "Almost," because we discuss only elasticity and fluid flow; we only 
touch on plasticity, and we neglect electric charges, magnetic coupling, and others. 
As stated in Chap. 1, a wide range of applications of phase-field modelling and 
simulations lies in this field. The models applied there, however, are complementary 
in concept with what we will discuss now: coupling of phase evolution with elastic 
distortion. 

It is reported that Armen Khachaturyan developed his branch of phase fields, 
the so-called “time-dependent Ginsburg-Landau theory,” to predict the equilibrium 
shape of a martensite needle in parent austenite. The problem relates to the old 
Eshelby problem of an inclusion in a matrix with finite transformation strain [9]. 
We may also include deviatoric strain if the crystal lattice changes from BCC to 
FCC. This was used to find the solution to the problem of the equilibrium shape of a 
crystal in its melt using the anisotropy of the interface energy in Chap. 3. Practically 
speaking, and this is often reported in the literature (to which I give no reference 
here), you cut out a piece of material from a crystal, transform it to a different phase 
with volumetric and deviatoric distortion with respect to the original crystal, and put 
it back. Continuity demands that both the parent and child crystal have to deform 
elastically to form a common body (see Fig. 7.1). 

In the phase-field approach, this is an easy exercise: you simply have to add the 
elastic bulk free energy of the individual phases to our previous models and decide 
on a good way to handle the elastic contribution within the interface. The elastic 
bulk free energy f* is (and all phase-field models agree with classical models 
from continuum mechanics on this point): 
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Phase transformation 
Uniaxial lattice expansion 


Inclusion and matrix under stress 


Fig. 7.1 Scheme of the 
elastic equilibrium between 
a transformed inclusion and 
the matrix 


1 
elast _ ^ 
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Y palei — CUM e — etl. (7.1) 
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where ej is the total strain in phase a, e*/J is the eigenstrain (or transformation 


free strain), and ci he is the elasticity matrix. We use the sum convention for double 
indices. » 

In general, e and C% kl are concentration- and temperature-dependent quanti- 
ties, which gives rise to intriguing effects of chemo-mechanical coupling (see [16] 
and references therein). We leave this for further reading. The ansatz (7.1) is a direct 
extension of the original multi-phase model for diffusive phase transformations, as 
the total elastic energy is a linear summation of the elastic energies of the individual 
phases weighted by the phase densities ø4. It is mostly assumed that the equilibrium 
is instantaneous. One derives the mechanical equilibrium equation, as we do with 
all relevant transport equations, as a functional derivative from the free-energy 
functional 


Davee F (1.2) 
de 

in vector notation, where the stress o and strain € are rank-two tensors. The difficulty 
now is to define these tensors within the interface, because the definition of the 
elastic free energy density (7.1) uses different strains €4. In all bulk phases, the 
stress is of course 0 = «C, but also the elasticity C is only defined for each 
phase. This is a classical problem in solid mechanics: to define an effective material 
from a mixture of materials, so-called mathematical homogenization. To be able to 
solve the mechanical equilibrium equation (7.2), one first has to define the effective 
mechanical properties of the interface, the effective elasticity matrix C, and also 
the effective transformation strain €*. The latter is commonly taken as a weighted 
average €* = »^, «€a. In the phase-field literature, the most-used homogenization 
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for the total strain tensor € is the so-called Voigt-Taylor model, which assumes 
homogeneous strain in all phases € = €4 = eg. Then, however, the stress has to 
be discontinuous in general. The opposite, the Reuss—Sachs model, takes the stress 
as continuous between the phases. Since the stress is the equivalent to a chemical 
potential in a mechanical system, one may take this assumption as an analogue to 
equal chemical potential in the interface (see also the discussion in [20]). 

A combination of both models, taking the strain as continuous in the tangential 
direction and the stress as continuous in the normal direction to the interface, is 
proposed in [8]. The most advanced model is the so-called rank-one convexification 
applying a jump condition (see [15, 17]). We leave this topic to further reading. Only 
one important consequence shall be detailed here: the elastic driving force for phase 
transformations Ag*!*' for a phase-transformation under mechanical load, either 
external or by internal stresses. For the Voigt-Taylor model, we find for the dual 
o—p interface (the general form is defined as a superposition between dual driving 
forces in the multi-phase-field approach): 


Agia = (= = a F 
B dha dog 


= (e — (Ee — €*)C + je — €*)? (Ca — Cg). (7.3) 


For Reuss-Sachs, we have: 


Ag = (e — e*)C(es — e5)C [e e$) ; (c =) Cle zi 
(7.4) 

We see in both cases, that the elastic driving force for phase transformations 
vanishes for a homogeneous material if e; = €? AND C, = Cg. Also in both cases 
we have two contributions: one related to the difference in the eigenstrains, the other 
related to the difference in the elasticity matrices. And of course, the elastic driving 
force vanishes if the total strain matches the eigenstrains, i.e., there is no stress in the 
interface region. Higher-order schemes considering a jump condition between stress 
and displacement at the interface [8, 15, 17] are more involved, but they follow the 
same principles. 

Another remark relates to the solution of the mechanical equilibrium equation 
(7.2). This can be achieved by any appropriate numerical approach using a finite- 
element scheme in real space. In phase-field models, the Fourier method is very 
popular because its solution for homogeneous systems, i.e., C — Cy — Cg, is com- 
putationally very cheap. If the material cannot be approximated by homogeneous 
elasticity, one first solves a homogeneous problem and then corrects the solution 
for the inhomogeneity in an iterative scheme. The difference between the elasticity 
coefficients is called "contrast." In the case of a low contrast, i.e., if the elasticity in 
all phases is similar, these schemes are quite effective. For high contrast, e.g., a pore 
filled with gas or liquid in a solid metal, one needs more elaborate schemes [12]. 
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One also has to consider that the Fourier transformation relies on periodic boundary 
conditions. In any case, however, one should not confuse the theoretical setup of 
the phase-field model with a special solution procedure, as important as this is for 
application. 

The final comment regards plasticity. There may be an outer load that exceeds 
the yield point of the material. Then, elastic strain is limited, and plastic strain sets 
in. This is associated with dislocation activity and the generation of another contri- 
bution to the free energy of the system: stored plastic energy related to dislocations. 
Such stored plastic energy may lead to phenomena such as “recrystallization.” This 
can be treated similarly to a phase transformation, and it is driven by the reduction 
of plastic energy. This is “easy” from a phase-field perspective, but it is difficult in 
terms of handling plastic energy at the continuum scale. It is even more difficult 
to design a model accounting for how moving interfaces interact with dislocations. 
Some applications are suggested for further reading. Additionally, internal stresses, 
which are caused by transformation strain between different phases, may exceed the 
yield point of the material, in particular at elevated temperatures. The stress state 
of the interface, and thus the driving force for phase transformation, is then limited 
by plastic relaxation. Therefore, plasticity will generate new driving forces on the 
one hand, and it will reduce driving forces on the other hand. This is a wide area 
for future research, but it is associated with special challenges regarding (i) model 
formulation and (11) numerical solution. 


Example—Martensitic Transformation 

A body-centered-tetragonal martensite structure is formed by rapidly quench- 
ing from the face-centered-cubic austenite phase. When a material undergoes 
a martensitic transformation, the crystalline unit cell undergoes a shape 
change. The final microstructure of martensite is strongly influenced by elastic 
energy generated during the change of the crystal structure. Accommodation 
of the elastic energy is often achieved by forming a multi-variant domain. 
In addition, plastic deformation occurs in the austenite matrix and the grow- 
ing martensite phase, and these are called self-accommodation and plastic 
accommodation, respectively [5, 23]. Thus, both elastic and plastic effects 
play important roles in martensitic transformations. 


The Simulation 

A 3D multi-phase-field model was used to model martensite microstructure 
formation in low-carbon steel. In this example (Fig. 7.2), an elasto-plastic 
approach is employed with a full set of 24 Kurdjumov-Sachs symmetry 
variants of martensite with real transformation strains [13]. A finite strain 


(continued) 
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Fig. 7.2 Lath martensite microstructures: (a) electron backscatter diffraction map of martensite 
in a sample containing 0.1 wt.% carbon [18]; and (b) simulated martensite microstructure in low- 
carbon steel consisting of 24 Kurdjumov-Sachs symmetry variants, indicated by color coding 


framework is applied in connection to a phenomenological crystal-plasticity 
model. Both methods are fully coupled within the OpenPhase library. 


7.2 Coupling to Fluid Flow 


In previous lectures, we have discussed the coupling of the phase field with 


temperature, solute, and elastic distortion. In all these cases, the driving force 


for a phase transformation Agag = (å — 55;) F is directly dependent on 


temperature, composition, or stress. We consider this coupling as “strong.” In 
the present subsection, we consider fluid flow, or melt flow when considering 
a solidification process. The coupling to phase evolution may be considered as 
“weak,” since the direct effect of flow on a phase transformation, the Clausius— 
Clapeyron effect, is generally very weak. The Clausius-Clapeyron effect describes 
the boiling temperature of a liquid dependent on the pressure in the system: the 
boiling temperature of water is significantly reduced on high mountains. The 
melting temperature of a metal, however, is hardly affected by the pressure of a 
streaming melt, and we therefore simply neglect it. There is, however, a strong effect 
of melt flow on dendrite morphologies in solidification. This is because on the one 
hand, melt flow significantly affects the transport of solute and heat in the melt; 
on the other hand, solidification significantly affects the viscosity of the material. 
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Fig. 7.3 Schematic of channel flow in which the width of the channel and the width of the interface 
are of the same order of magnitude 


A solid dendrite can be seen as a rigid body. If it is attached to the mold of the 
casting, it forms a rigid barrier for flow; if it is transported with the melt as an 
equiaxed crystal, it will affect the effective viscosity of the two-phase system (solid 
and melt). If the solid fraction is low, a metallic melt has a viscosity like that of 
water; if the solid fraction exceeds 30%, it will behave like honey with precipitated 
sugar crystals; above 5096, any melt flow will stop. 

This mutual interaction leads to intriguing phenomena, which we will leave 
for further reading. Here, we will elaborate on one effect that is crucial for the 
interaction of a fluid with a solid when the interface between liquid and solid is 
diffuse at a mesoscopic scale. At the microscopic scale, we definitely accept a 
so-called “no-slip” condition: the flow velocity decays to 0 monotonously in the 
direction normal to the interface. Figure 7.3 schematically shows the condition of 
a channel flow in which the width of the channel and the width of the interface 
are of the same order of magnitude. One side of the channel is treated as a “wall,” 
i.e., a sharp interface with a well-defined no-slip condition; the other side is treated 
as diffuse in the context of phase-field theory. How do we realize an analogue to 
the no-slip condition within a diffuse boundary? As in the first part of this lecture, 
we treat the interface as an effective material, and in the present case as a porous 
medium: the flow will penetrate the interface, but the permeability of the interface 
will be a function of the phase field. 

The corresponding fluid-flow equation for the fluid velocity u in a mixed 
domain—solid and liquid—is readily written down (for simplicity, without moving 
solids, i.e., the solid velocity is set to 0; @ = riquia, density p set to 1): 


ER E E å 
ay oH + Vouu = SVP + V (vVøu) — h* Xj. (7.5) 
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The important part here is the last term: friction within the diffuse interface Xj. 
This has been worked out for the planar interface with a friction coefficient h*.! 
Its numerical value is given as h* = 2.757 [6]. Other models defining the no-slip 
condition as a function of the phase field have been investigated [3, 4, 21] and may be 
investigated in the future. A consistent investigation for curved interfaces (concave 
or convex) is still missing. In any case, we will require the “sharp-interface solution” 
to be met outside of the interface, i.e., that we meet the sharp-interface model in the 
bulk. In the case of the channel flow, the sharp-interface model predicts a parabolic 
velocity profile in the channel: Hagen-Poiseuille law. The maximum velocity is 
a sensitive indicator of whether the diffuse interface correctly emulates the sharp 
interface; it is a kind of “thin-interface limit,’ i.e., we match a sharp-interface 
solution in the bulk liquid but deviate systematically within the thin interface. 


7.3 Exercises 


Exercise 
Derive the expressions for the elastic driving forces (7.3) and (7.4) for the 
Voigt-Taylor and Reuss-Sachs limit, respectively. 


Example: Dendritic Solidification Interacting with Shear Flow of the 
Melt 

The setup (Fig. 7.4) represents a thin Mg-Al melt channel between two rigid 
walls. We see nucleation and growth of a-Mg dendrites. At the same time, 
shear flow is introduced to the melt via the motion of the rigid walls. Full 
integration of the inertial and friction forces acting on the solid dendrites and 
the melt results in the dendrites moving with the melt. In this example, the 
fluid-flow problem is solved using the lattice Boltzmann method, while the 
system morphology and its evolution is described by the phase-field method. 
Both methods are fully coupled within the OpenPhase library, allowing the 
study of arbitrarily complex geometries. See [14, 22]. 


! h* is named the “Hermann-Joseph” constant. Hermann-Joseph Diepers was one of my early 
collaborators. He passed away from cancer before finishing his PhD. His memory lasts. 
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Fig. 7.4 Flow simulation of Mg-Al alloy solidification 


Further Reading 


* Chemo-mechanical coupling [10, 16]. 

* Hadamard jump [15, 17]. 

* Recrystallization and rafting under high-temperature creep [1, 2, 7, 11]. 
* Dendritic growth with buoyancy [19, 22]. 
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Chapter 8 A 
Quantum Phase Field FEEN 


8.1 Introverted Picture of Mass and Space 


How to Start? There are so many aspects coming together to form a consistent the- 
ory, culminating in phase-field theory. From the physical point of view, this includes 
thermodynamics, wave mechanics, interface-driven phenomena and capillarity, and 
the kinetics of phase transformations. From the practical point of view, we have 
the numerics of the solutions of partial differential equations, and multi-physics 
problems. The main applications of phase-field theory to date have been engineering 
problems. From the mathematical point of view, we derive a set of partial differential 
equations using variational principles from a governing functional that has a 
gradient contribution acting on the field. In quantum-mechanical language, this 
gradient is a momentum operator, and the field is a quantum field with a discrete 
spectrum. If this sounds strange, then be reassured that, as before, we will go through 
this step by step. 


Let us start from the conceptual point of view. We discuss a free-boundary 
problem. The free boundary is the inner boundary between two or more domains, 
which are characterized by phase fields. In classical phase-field theory, we call the 
boundary an interface. The boundary is termed “free” because it is not fixed by 
some external condition; it can evolve in time (see Chap. 1). This inner boundary 
is different from outer boundaries, i.e., the boundaries of the calculation domain, 
which are fixed. The inner boundary between phases evolves in space and time. 
This evolution is closely coupled to transport phenomena in the domains: transport 
of temperature, solute, and momentum, as well as displacement. These phenomena 
have been discussed in the previous lectures. 

The object of interest in phase-field theory is the inner boundary. This is not 
something that we specify; it is something that will come out: it “emerges.” It shall 
be independent of the settings of the outer boundary in any system; if there is any 
outer boundary at all. 
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Fig. 8.1 The billiard table as a representation of the extroverted view of space in traditional 
physics: objects are placed in space. The network of introverted spaces (yellow lines) connecting 
the balls reverts the view to a closed network system without outer boundaries (© www.winsport. 
de with permission) 


Boundaries bound “spaces”. We may distinguish between “introverted” and 
“extroverted” spaces. The latter means the space around an object. The extroverted 
space is a kind of background, existing with or without the object. This is the 
traditional view of space. Objects are placed in space. As an illustration see Fig. 8.1. 
The green of the billiard table represents space. The objects, billiard balls, are placed 
in this space. They interact according to Newton’s laws, momentum and energy 
conservation. In a real billiard game the player uses the cushions to mirror space. 
We may, hypothetically, push the cushions to infinity, then our play ground will 
be infinitely large, but balls would never come back. We may postulate periodic 
or “Moebius” boundary conditions. But all this becomes irrelevant if we revert 
the picture to an introverted space. Introverted space lives inside the objects: it 
is bounded by the objects. These introverted spaces are sketched by the yellow 
lines in Fig. 8.1. The relative position of all balls is uniquely determined by the 
network of inner spaces. An absolute outer space is not needed to characterize their 
interrelation. 

In the physical world, which is composed of particles and space, the only objects 
that could serve as a bound of an introverted space are elementary particles. So, if 
particles (zero-dimensional objects) are the boundaries, space is the 1-dimensional 
line-distance between particles. There is no other choice. 

If there are many particles, there will be many spaces. In the quantum-phase- 
field concept, each space is associated with a phase field $,. From the number 
of elementary particles within the observable universe N, a maximum of N 
combinatorial 2 spaces is a number of the order of 10!??. Space will also be taken as 
“substantial,” i.e., it is not just “empty”: it has an energy, a negative energy related to 
vacuum fluctuations, as shown in Sect. 8.3. The 3D “space of cognition,” the space 
of our daily experience, will be reconstructed later, in Sect. 8.4. We will proceed like 
this: Define a system of quantum phase fields, closed in itself, without having outer 
boundaries. 
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8.2 Formal Definition of Quantum Phase Fields 


We define (and you have all the background from the previous lectures) a set of a 
large number of phases $4. The phases form a system that is closed in itself, forming 
a "universe." It is a system without outer boundaries. We do this in analogy to the 
multi-phase-field system from Chap. 6: 


Va =1. (8.1) 


Here, $4 is simply termed the “phase” since up to now no “space” has been 
defined. Each phase ø4 is an element of the universe that is different from all other 
elements $g. Therefore, an interface between these phases is formed: a fermionic 
particle, as we will see. The phases incorporate the basic elements of physical 
matter, “mass” and “space.” We will associate them with the conserved quantity: 
energy H. The previous statement that "the phases incorporate mass and space" is 
important: the phases (our phase fields) $4 are not defined in space, as in classical 
or quantum field theories, but rather they define space! Massive elementary particles 
are not “placed into space,” they “connect spaces.” We will derive this step by step. 

The energy functional of the system in quantum mechanics Å is an operator 
(as indicated by the hat). This is acting on the wave function of the system |w). In 
our case, the system—though this sounds quite highbrow—is the universe; |w) is 
the wave function of the universe. The only ingredient of the theory, up to now, is 
energy H = (w| Å |w) as the expectation value of Å. Energy is conserved (first 
law of thermodynamics).! Furthermore, we will postulate that H = 0,i.e., there is 
no net energy: all energetic states, positive and negative, have to sum up to 0. The 
argument for this is simple: there is no evidence regarding where a finite energy 
should come from (see also the Wheeler-DeWitt theory [7]). 

We will allow changes d Å, that the zero-energy state, the state of "nothing," 
separates into positive and negative energetic states, the state of “something.” This 
may be related to the Big Bang as the origin of our universe, if you will. Changes 
are related to a non-conserved quantity, entropy: the second law of thermodynamics. 
We will expand H, in the changes d Å with respect to the phases øy. Å thereby is 
itself a function of all fields (ø,), H=H ({¢g}), and, as a reminder, all fields are 
connected by the sum constraint (8.1): 


å 1 0H 
a=y f dør (8.2) 


! Note that the constraint of energy conservation is in contradiction to general relativity [8], where 
energy is not conserved! We will accept this. 
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The integral runs over the definition range of the phase fields from 0 to 1, meaning 
yes or no, existing or not-existing, and we allow the fields to vary between these 
bounds, i.e., that they are diffuse, as is usual in phase-field theory. 

The phases will be used to indicate different states of energy. This means that they 
will be used to derive relations between the different states of energy. We will allow 
a BIG number of fields connecting a BIG number of different states of energy, both 
positive and negative. An important factor in this regard is that positive and negative 
states are not just “mirrored” states like matter and anti-matter in traditional physics: 
they are topologically different, as we will discuss in closer detail when we discuss 
the ordering scheme of space and mass in Sect. 8.4. For now we simply state, that 
there are negative energy states Ey associated to a phase o, and positive states Ugg 
associated to junctions between phases o and £. 

There is no fundamental space. We repeat the argument: where should it come 
from? 

“Space” is introduced as a distance Qg, defined by the inverse of the negative 
state of energy Ey < 0 of the field ¢y: 


MER ur (8.3) 
in RE,” ) 


where h is Planck's constant, c is the speed of light, and à is a positive and dimen- 
sionless coupling parameter to be determined. It will be proven self-consistently in 
Sect. 8.3 that this definition leads to a 1D line coordinate sg specific to each phase, 
du = $a (sa). Having this line coordinate we can treat Øy (sy) as a classical phase 
field in 1D. 

There is no proof that this is the only possible construction to define our universe, 
no rigorous argument against concepts of “parallel universes” or “‘multiverses,” or 
10- or 21-dimensional spaces with strings and branes embedded, as postulated by 
several researchers (see, e.g., [14] for 10 dimensions and [18] for 21 dimensions). 
You will find almost every number of dimensions for space-time postulated from 4 
to 21, maybe even higher, since the original work of Kaluza and Klein [15, 16]. 
A whole discipline, called “mathematical physics,” is searching for a “unitarian 
description of all physical phenomena” in multidimensional manifolds with appro- 
priate geometrical measures. 

We will approach this differently. The present theory introduces time and space 
as auxiliary coordinates: the coordinates are not “fundamental.” It is shown that 
the theory, which is based on a set of principal statements, is consistent with these 
statements when formulated in these auxiliary coordinates. It is “self-consistent,” 
and it resembles our ee of the physical world. 


Substituting døy = Tte dsa and introducing the forces ha = gå Å vields 


N 
> © Oda Pee TR 
a 2. pa a -X dsah({$p}). (8.4) 
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Fig. 8.2 Two doublons in a ó 
periodic setting. Each 

doublon is formed by a l 
right-moving and a 

left-moving soliton. The 

velocity is proportional to the 
energy difference between 
adjacent doublons 


Space emerges, i.e., it is created by variations in the phases d$: compare the 
"diffuseness of phase field.” We relate the line coordinate s, to the distance Qw, 
defined by the inverse of the energy quantum Eg (see (8.3)), by the integral 


f. hc 
didis enc. (8.5) 
FT aPa a 48E, 


This is easily be seen from the solution of $5 (Sq), Fig. 8.2. Phases and space are 
complementary objects, defined by phase fields. In the present concept, the phase 
field generates space instead of living in a space. But we need to swallow a toad: 
space sy is a 1D line coordinate. There is no evidence of a 3D space (besides our 
daily cognition, which might be in error!). 

To construct the energy operator, or Hamiltonian H, we use the standard form of 
the Ginzburg—Landau functional in 2D Minkowski notation:? 


. “rte aun{f(a.\? 1/83 P r 
á-y[ ds 72 ED z2 (4.) + r Pad] 


a=1 


(8.6) 


You will recognize the similarity with the free-energy model with a double-obstacle 
potential from Chap. 2. In addition to the gradient contribution in space, we include a 
gradient contribution in time. The time derivative accounts for dissipation, i.e., that 
in contrast to gradients in space, which have a positive energy penalty, gradients 
in time are favored energetically: the system shall evolve and not stagnate. We 
may treat this philosophically or simply state that this ansatz ensures relativistic 
invariance due to Lorentz contraction of the interface width (see [20] for details). 


? We have removed the index o from space and time coordinates for readability. 
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U is a positive energy quantum to be associated with the positive rest-mass of 
elementary particles, the equivalent of the surface energy in the traditional phase- 
field model? The gradient contributions of the Hamiltonian, Eq. (8.6), x and ig, 
shall be understood as operators acting on a quantum mechanical wave function 
|w). This allows to evaluate the expectation value of the energy for an actual state 
of the system. The wave function in “quasi-static approximation” will be explicitly 


constructed below in Sect. 8.3. The phase-field equation is written down: 


a [o H 27 
tafe = — 5G, : (w|H|w). (8.7) 


Since gradients of time are considered, we have to integrate over time to have 
a well-defined functional derivative, as discussed in Chap. 2. This equation has 
two parts: (i) a non-linear wave equation for the phase field; and (ii) a linear 
Schródinger-type equation for quantum-mechanical excitations. This procedure is 
not new; it can be traced back to the so-called de Broglie-Bohm double-solution 
program [1, 2, 5, 6]. It is fully accepted as an alternative interpretation of quantum 
mechanics compared to the prevailing so-called “Copenhagen interpretation.” For 
more explanation, see [17]. 

Now, we separate the expectation value of the energy functional (8.6) into three 
different contributions. These are distinguished by whether the differential operators 
E and 2 are applied to the wave function |w) or the field ø4. 

Applying the differential operators to the phase-field components and using the 
normalization of the wave function (w|w) = 1 yields the force ug [I] related to the 
gradient of the fields o: 


|. 4Un| (ØRN 1 (ØRN 7? 
Ua = 72 ( as ) c ( Jt ) + 72 (Pa tol] (8.8) 


This contribution relates to the interface energy in the conventional 3 dimensional 
phase field application. In 1 dimensions it is a force quantum, related to gradients of 
the phase field or to values of the phase-field 0 < øy < 1. 

The mixed contribution, when one of the operators 2 and È is applied to the 
field ¢ and one to the wave function |w), describes the correlation between the field 
and the wave function. It shall be set to 0 in the quasi-static limit. In this limit, we 
keep the field static for the evaluation of the quantum-mechanical force. Then, we 
take this force for the determination of the time evolution of the field. A coupled 
solution has not been worked out to date: 


Eam 3 1 3øa Eu 


72 ðs aTa ga 


0 = (1— 249) (8.9) 


? Note that the unit here is Joule, J, instead of 
living in a 1D space. 


ES for the surface energy in 3D, because we are 
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It is shown in [20] that this Eq. (8.9) is consistent with Newton's second law of 
a 


acceleration. Finally, we apply the momentum operators 3- and 73. to the wave 


function |w), which yields the force ej [I]: 


32 13 


4Un 
E as? 3t n 


€g = 


a. (w| (8.10) 


72 

This contribution applies to the bulk energy of the phase field ø4 = 1. We will 
explicitly evaluate this after the structure of the solutions of the fields is discussed. 
For a steadily moving field with velocity v, one transforms the phase-field equation 


into the moving frame traveling with this velocity è = vå. Inserting Eqs. (8.8)- 
(8.10) into (8.7), we find: 


0 ô Too 
T— Øy = ——— dt(w|H 
Fete == 5 C 
825, v? n? 1 
=U L 352 (1 5) + i (o. 3) + mg, Ae, (8.11) 
where: Ae = eg — eg is the difference in the volume force (the equivalent of 


the Gibbs free-energy difference) between two fields, Eq. (8.10); and m, is the 
appropriate coupling function. As a first result, we see that the interface width 


n contracts with velocity: ny = ,/1— 5, Lorentz contraction. It collapses for 
v — c, and velocities v > c are forbidden. 

The structure of the field is well known to us from Chap. 2: the minimum solution 
of Eq. (8.11) is the “soliton.” A periodic solution for two fields is displayed in 
Fig. 8.2. One field bounded by two solitonian waves, one right-moving and one left- 
moving, is called a “doublon.” This is the basic element of physical space øy = 1, 
bounded by the junction to the other field, which will be interpreted as an elementary 
particle later. 


8.3 Volume Energy of One Doublon 


From the doublon solution in Fig. 8.2, we can see that the field forms a 1D box with 
fixed walls and size €2, for the field œ. This is a standard exercise for a quantum 
problem: the particle in a box potential. The difference here is that the particle is 
massless, i.e., the dispersion relation is linear in momentum p instead of quadratic 
in p as for the massive particles. According to o Casimir [3], we have to compare 
quantum fluctuations in the box with discrete spectrum p and frequency wp = a 
to a continuous spectrum. This yields the negative energy Eq of the field a: 


h o0 h 
Ecc. ye-f pdp | = -å £ (8.12) 
1 
p= 
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where à is a positive, dimensionless coupling coefficient. I have used the Euler- 
Maclaurin formula in the limit e — 0 after renormalization p > pe *P. Since all 
parameters are positive, we see that "space" is accounted for by negative energy, 
scaling inversely proportional to the size of the doublon, which proves (8.3) for 
self-consistency. This energy scales like the energy of the gravitational field in 
Newtonian mechanics. Therefore, the force, as the derivative of the energy (8.12) 
with respect to space, can be associated with a gravitational attraction between 
the junctions between different doublons. The junctions can thus be interpreted 
as massive objects: elementary particles. They are associated with positive energy 
U [Eq. (8.8)] corresponding to the interface energy of a traditional phase field, 
while the bulk energy of the doublon is negative. In contrast to Newtonian 
mechanics and in agreement with general relativity (see "gravitational waves" [9]), 
the attraction is a wave phenomenon with time-dependent action. The coefficient a 
can be determined from the measured gravitational constant on Earth. Applying 
the theory to all masses of the visible universe gives a prediction for repulsive 
gravitation at ultra-long distances and an explanation of the expansion of the visible 
universe [20, 21]. 


8.4 Multidimensional Interpretation 


As stated at the beginning, the present concept has no fundamental space. The 
distance $2, is intrinsic to one individual doublon o. Different doublons are 
connected in junctions due to the sum constraint (8.1): a multiple junction in the 
multi-phase-field terminology defines an elementary particle in the quantum-phase- 
field concept. The position of one particle related to an individual component of the 
field is determined by the steep gradient Og *. In a multiple junction, where many 
fields intersect, the junction couples many different directions related to the different 
doublons. We shall consider the junction (or particle) as a “very small volume,” in 
physics terminology, this is a zero-dimensional object. Since half-sided solitons are 
spinor-type objects belonging to the 3D SU(2) symmetry group, we may order the 
incoming and outgoing doublons of a junction in 3D Euclidean space. We will call 
this space the "space of cognition," since our cognition orders all physical objects 
in 3D space (Fig. 8.3 sketches this picture). Individual doublons form a “doublon 
network.” Each doublon is expanded along a 1D line coordinate and bound by two 
end points, described by gradients of the field, right- and left-moving solitons. Due 
to the constraint (6.6), the coordinates of different fields have to be synchronized 
within the junctions of small but finite size n. The constraint (6.6) also dictates that 
there are no “loose ends"; the body is closed, in itself forming a “universe.” 
Finishing this lecture let me remind you: In the language of network theory, the 
junctions of a network are called vertices (in Latin the plural of vertex, a single 
junction or knot). The connections between vertices, are called “edges”. Such a 
network can be embedded into a 3-dimensional vector space, or higher (but not 
lower). The doublon network, as we call it in the context of this lecture, does not 
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define a vector space by itself! This means that empty spaces inside the meshes 
of the network are not accessible, as they would in a vector space. They have no 
physical reality. Only points on the doublons can be attributed by a local coordinate, 
related to the line distance to a junction or vertex. This is important, since vertices, 
which are not connected by edges, are "invisible:” there is no physical space 
between them where light, information, action etc. could be transmitted through. 
It will be a future task to embed the doublon network into a four-dimensional space- 
time where the doublons form geodesic trajectories between particles. 

A last comment relates to the topology of positive and negative states of energy. 
The positive states U relate to particles, junctions between the doublons: each 
particle connects all incoming and outgoing doublons. In contrast, one doublon 
connects exactly two particles. This defines the network structure, see also above. 
Furthermore, positive and negative energy states cannot simply annihilate without 
changing the topology of the whole system: once separated into MANY positive 
and negative states, the system evolves irreversibly. 


8.5 Symmetry Breaking in Condensed Matter and 
Elementary Particle Physics 


"Symmetry breaking" is a fundamental concept in materials physics, see the consid- 
eration about “phases” and the “order parameter" in Chap. |. The term “symmetry” 
here relates to solid materials’ crystal symmetry, or the lack of symmetry in the 
amorphous state, liquid or gas. A particular case is magnetic phase transformation 
from the ferromagnetic to the paramagnetic phase. In the ferromagnetic state we 
see spontaneous magnetization of the magnetization +M , which is a vector in 3- 
dimensional space and can have + or — direction. This is the “symmetry broken 
state”. In contrast, the paramagnetic phase, which is the stable phase above the 
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Curie temperature To. It is called the “symmetric state" with vanishing spontaneous 
magnetization M =0. 

In the context of elementary particle physics, the concept of "symmetry break- 
ing" was first introduced by Goldstone [11], referring to the theory of superconduc- 
tivity, i.e. a phase transformation in condensed matter. Very soon, the publication 
"Broken Symmetry" by Goldstone et al. [12] appeared, which established the 
concept of symmetry breaking in elementary particle physics, so-called Goldstone 
modes. Only two years later, the work of Higgs [13] and Englert and Brout [10] 
appeared, which was awarded the Nobel prize in 2013 after the so-called “Higgs 
boson" had been confirmed experimentally [4]. In light of the present concept, the 
Higgs field, which couples to fermions to give them their mass, can be identified 
with the phase field as we call it in materials science: An order parameter field in 
space and time with broken symmetry. The underlying mathematical formalism, 
Lagrange and Hamiltonian functions are almost identical. The new contribution of 
the “quantum phase field" is the explicit solution of the minimum energy solution 
of this field in 1 dimension; the doublon, as detailed. Doublons then are connected 
to form a doublon network in 3+1 dimensional space-time by the sum constraint of 
a multi phase-field, Eq. (8.1). 


yc? =1. 


In conclusion: phase-field theory has two different roots (see Chap. 1): one 
in thermodynamics and one in wave mechanics. Thermodynamical principles of 
energy conservation and entropy maximization set the stage. The wave-mechanical 
picture of phase fields rests on the theory of solitons. Quantum-phase-field theory 
has roots in de Broglie and Bohm's interpretation of matter as a wave phenomenon; 
it combines the soliton solution of a non-linear wave equation, the phase-field 
equation, with a Schródinger-type linear wave equation in a finite domain. The 
non-linear wave equation is thus derived from the thermodynamic theory of phase 
transformation: Ginsburg-Landau theory. We may use other potentials than the 
double well or double obstacle, as explained earlier. We may use higher gradients in 
o. The general structure, however, is consistent, having three different contributions, 
as in every phase-field model: gradient, potential, and bulk (see Chap. 2). 

Where is the bulk part in the Hamiltonian (8.6)? The bulk part arises from 
applying the gradient operator to the wave function! It is defined from the solution 
of Schroedinger equation in the box formed by the doublon. In the concept, 
and in “reality” (I think), there is no given "space," no bulk in which objects, 
"elementary particles," are placed. The elementary particles, which in the phase- 
field interpretation are junctions or knots between doublons, are defined by gradients 
of the phase fields. They are the endpoints of the doublons and form the bounds of 
space. Space and mass are two sides of the same coin. This all may seem "quite 
philosophical" to the student, or simply “confusing.” But it should not discourage 
you from doing the exercise below [calculating the energy of “space” from (8.10)]. 
We conclude with the following statements: 
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Energy is substantial and conserved. 

There is no fundamental space, no fundamental time. 

Space, a one-dimensional distance, emerges from variation of energy. 

Mass is attributed by positive energy, space by negative energy. 

Two half-sided solitons form the doublon, which is the primitive object of the 
concept, defining mass and space. 

The doublon belongs to the 3D S(U2) group of spinors. Doublons span the 3D 
space of our cognition. 

The concept is relativistically invariant and makes predictions at the scale of the 
universe (see "further reading"). 


8.6 Exercises 


Exercise 
Derive the “energy of space" (8.12): E, = —å så defined by vacuum 
fluctuation in the 1D box formed by a doublon of size $25. 


Further reading 

* Hydrodynamic basis of quantum mechanics and its relation to Phase 
Field: [19] 

* Expansion of the universe by repulsive gravitational action on ultra-long 
distances: [20, 21]. 


References 


1. 


2. 


3: 


5. 


6. 


D. Bohm, A suggested interpretation of the quantum theory in terms of “hidden” variables. I. 
Phys. Rev. 85, 166—179 (1952) 

D. Bohm, A suggested interpretation of the quantum theory in terms of “hidden” variables. II. 
Phys. Rev. 85, 180-193 (1952) 

H. Casimir, On the attraction between two perfectly conducting plates. Proc. Koninklijke 
Nederlandse Akademie van Wetenschappen B51, 793—795 (1948) 


. S. Chatrchyan, et al., Observation of a new boson at a mass of 125 GeV with the CMS 


experiment at the LHC. Phys. Lett. B 716(1), 30-61 (2012). ISSN: 0370-2693. https://www. 
sciencedirect.com/science/article/pii/S0370269312008581 

L. de Broglie, Nonlinear wave mechanics, in Trans, ed. by A.J. Knodel (Elsevier, Amsterdam, 
1960) 

L. de Broglie, L'interpretation de la mechanique ondulatoire par la theorie de la double 
solution. Proc. Int. School Phys. Enrico Fermi 49, 346—367 (1971) 


90 


8 Quantum Phase Field 


. D.S. DeWitt, Quantum theory of gravity. I. The canonical theory. Phys. Rev. 160, 1113-1148 


(1967) 


. A. Einstein, Die Grundlage der allgemeinen Relativitåtstheorie. Ann. Phys. 49, 769—822 (1916) 
. A. Einstein, Über Gravitationswellen. Sitzungsberichte der Kóniglich Preussischen Akademie 


der Wissenschaften Berlin 154—167 (1918) 


. F. Englert, R. Brout, Broken symmetry and the mass of gauge vector mesons. Phys. Rev. Lett. 


13, 321—323 (1964). https://link.aps.org/doi/10.1103/PhysRevLett. 13.321 


. J. Goldstone, Field theories with « superconductor » solutions. Il Nuovo Cimento 19(1), 154— 


164 (1961). https://doi.org/10.1007/BF02812722 


. J. Goldstone, A. Salam, S. Weinberg, Broken symmetries. Phys. Rev. 127, 965-970 (1962). 


https://link.aps.org/doi/10.1103/PhysRev.127.965 


. PW. Higgs, Broken symmetries and the masses of gauge Bosons. Phys. Rev. Lett. 13, 508-509 


(1964). https://link.aps.org/doi/10.1103/PhysRevLett.13.508 


. L. Hughston, W. Shaw, Classical strings in ten dimensions. Proc. R. Soc. A Math. Phys. Eng. 


Sci. 414, 423-431 (1987). https://doi.org/10.1098/rspa.1987.0152 


. T. Kaluza, Zum Unitåtsproblem der Physik. Sitzungsberichte der Königlich Preußischen 


Akademie der Wissenschaften 966—969 (1921) 


. O. Klein, Quantentheorie und fünfdimensionale Relativitütstheorie. Z. Phys. 37, 895—906 


(1926). https://doi.org/10.1007/BF0139748 1 


. J. Kundin, I. Steinbach, Quantum-phase-field: from the Broglie-Bohm double-solution pro- 


gram to doublon networks. Z. Naturforschung 75(2a), 155—170 (2020). https://doi.org/10.1515/ 
zna-2019-0343 


. M. Loev, Origin of everything and the 21 dimensions of the universe, in APS March Meeting 


Abstracts. APS Meeting Abstracts, S1.104 (2009) 


. R. Mauri, A non-local phase field model of bohm's quantum potential. Found. Phys. 52, 52-58 


(2001). https://doi.org/10.1007/s10701-021-00454-9 


. I. Steinbach, Quantum-phase-field concept of matter: emergent gravity in the dynamic 


universe. Z. Naturforschung A 72(1), (2017). https://doi.org/10.1515/zna-2016-0270 


. I. Steinbach, J. Kundin, F. Varnik, Self similarity of the expanding universe as understood by 


quantum-phase-fields. arXiv: 2002.12848 [physics.gen-ph] (2020) 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, 
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate 
credit to the original author(s) and the source, provide a link to the Creative Commons license and 
indicate if changes were made. 


The images or other third party material in this chapter are included in the chapter's Creative 


Commons license, unless indicated otherwise in a credit line to the material. If material is not 
included in the chapter's Creative Commons license and your intended use is not permitted by 
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from 
the copyright holder. 


Part II 
OpenPhase 


Chapter 9 A) 
Tutorial 1: OpenPhase TCA 


9.1 Introduction to the OpenPhase Software Package 


OpenPhase is an open-source C++ software project started in 2008 at The Inter- 
disciplinary Centre for Advanced Materials Simulation, Ruhr-Universität Bochum, 
Germany. OpenPhase is dedicated to investigating microstructure evolution in mate- 
rials undergoing first-order phase transformations, capillarity-driven coarsening 
processes, or a combination of both. Also you may embed the simulation into 
a macroscopic process simulation which determines boundary conditions on the 
micro domain, heat fluxed, mechanical loading or others. 

OpenPhase is based on the multi-phase-field model as developed by the lecturer 
and his co-workers. However, the software is not restricted to “Steinbach models.” 
Other models, such as the Khachaturyan scheme of interface homogenization (see 
Chap. 7), are already implemented to be applied to problems where these models are 
of advantage, either from a theoretical point of view or simply from a practical point 
of view for an effective numerical solution. It is straightforward to implement your 
model in the software. We also strongly encourage scientists all around the world to 
provide their own special solutions. You are welcome to contribute new modules to 
the software project! 

The commercial version, OP Studio, is meant for users from industry, as well 
as for users from academia who do not have sufficient experience in numerical 
simulations to handle the academic code without support. It provides a graphical 
user interface and additional features like the coupling to thermodynamic databases 
for multi-component diffusion, as well as finite-strain elasticity and advanced 
plasticity models. Furthermore, custom solutions for specific user problems can be 
offered. 

The first section of this chapter is dedicated to providing an overview of the 
OpenPhase structure, including installation, compilation, and executing your first 
example. In the second section, we explore several representative examples that 
demonstrate the use of the software and its effectiveness. 
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9.1.1 Features 


e Open-source library. 

* Simplified microstructure creation. 

* |D,2D and 3D simulations are possible. 

* Multiple chemical components, phases, and grains are possible in a multi-phase- 
field simulation. 

e OpenMP and MPI parallelization. 

* Tools to extract various statistics data are provided. 


9.) Download 


OpenPhase is licensed under the GNU General Public License version 3 
(GPLv3) for all open-source applications. We invite you to visit our website at 
https://openphase.rub.de to obtain your OpenPhase version. 

Paid support and commercial licenses for the OpenPhase library can be obtained 
from the OpenPhase Solutions GmbH. Visit: https://openphase-solutions.com. 
OpenPhase Solutions GmbH also offers an end-user-oriented standalone phase- 
field application with extended functionality. 


9.3 Installation 


To install OpenPhase follow the installation instructions provided here. Make sure 
to manually install the required packages to ensure a successful installation. 


9.3.1 System Requirements 


OpenPhase works on Unix-like OSs. Therefore, if you would like to use Windows, 
we advise you to install the Windows Subsystem for Linux (following Microsoft’s 
instructions: — https://docs.microsoft.com/en-us/windows/wsl/install#manual-inst- 
allation-steps) or to install any Linux distribution as a virtual machine on Windows. 


Minimum System Requirements 


* GNU GCC C++17 compliant compiler (GCC @ 9.0.0 or greater). 
* The Fastest Fourier Transform in the West (FFTW) package 
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9.3.2 Installing the Prerequisites 


The following software is required before installing OpenPhase: 


* The GCC compiler (g++) 
e FFTW library. (www.fftw.org) 
The GCC Compiler 


1. Open the terminal (Linux: Ctrl+Alt+T) 


2. Update your packages using: 


$ sudo apt-get update 


3. Install the development package: 


$ sudo apt install build —essential 


4. Check g++ compiler version: 


$ gcc ——version 


FFTW 


1. If you are using the apt package manager, Installing fftw-dev package is as easy 
as running the following command on terminal: 


* Non MPI users 
$ sudo apt-get install fftw—dev libfftw3—3 libfftw3—dev 


* MPI users 
$ sudo apt-get install libfftw3—mpi—dev 


2. Or you can use the source code instead. Download the newest FFTW package 
using https://www.fftw.org/download.html. 


* Now navigate to the Downloads folder and extract using: 


$ tar xvzf FileName.tar.gz 


* Navigate to the extracted folder. The installation can be as simple as: 


$ ./configure 
$ make 
$ make install 


* The FFTW MPI source code can be compiled by following: https://www.fftw. 
org/fftw3 doc/FFTW-MPI-Installation.html. 
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9.3.3 OpenPhase (Compilation) 


After downloading the source code of OpenPhase as described in Sect. 9.2, extract 
the source file and navigate to the extracted folder. 
Now, compile the source code as follows: 


* Non MPI users 


$ make 


* MPI users 
$ make SETTINGSzmpi- parallel 


This will take a few minutes, depending on your system. 
If you got a ^Compilation done" message on your terminal, you are ready to run 
your first example. 


9.3.4 Update OpenPhase 


OpenPhase employs standard versioning, although it is continuously under devel- 
opment with frequent updates. We recommend keeping OpenPhase up to date 
as you use it to develop your application(s); regular upgrades are advised. New 
and archived versions of the OpenPhase library can be found on our website at 
https://openphase.rub.de/download.html. 


9.3.5 OpenPhase Application 

Structure of an OpenPhase Application 

Each application example directory should have at least the following files (Manda- 
tory): 


* ProjectInput.opi 
* Main.cpp 
* Makefile 


It may include additional files (Optional): 


* Auxiliary input files (e.g. EBSD, Orientations, Geometry, etc.) 
* Post Processing files (e.g. Gnuplot, python, etc.) 


The structure and parameters for the input file, which has been given the name 
"ProjectInput.opi", will be discussed in the next section. 


9.4 The Input File Structure 97 


Compile Your Application 


Compile the application source code as follows: 


* Non MPI users 


$ make 


* MPI users 


$ make SETTINGSzmpi- parallel 


If the application is functioning properly, the output won't indicate errors. 
Depending on your compiler version, you may get warnings that can usually be 
disregarded. Be sure to recompile and test your application after making changes to 
the source code or updating OpenPhase library. 


Run Your Application 


* Non MPI users 


$ /ApplicationName 


* MPI users 


$ mpirun —np 4 ApplicationName 


This command tells the computer to run the executable program "Application- 
Name" on 4 CPUs. 

Note: This mpirun command may be useful for novices with simple tasks, but 
it's not the way to go when you have a large, lengthy application! Additionally, 
outcomes can vary from machine to machine. When you request 4 processors, for 
example, 


* [n the best-case scenario, mpirun detects three additional PCs similar to the one 
you're using, copies your software to them, and runs it on all four. 

* [naless ideal scenario, your current system has four processors, this memory is 
shared among them, and your program runs; 

* In the worst-case scenario, there are less than four processors available, therefore 
one of them will “timeshare” and alternately run two or more of your process. 


9.4 The Input File Structure 


OpenPhase is a command-line, non-interactive tool; it anticipates that the simulation 
input parameters will be provided in the form of an input file. The user must ensure 


98 9 Tutorial 1: OpenPhase 


that the input parameters file is placed in the same folder as the simulation example 
before proceeding. This file is typically referred to as ProjectInput file. The 
OpenPhase software will read the input file and begin the simulation immediately if 
the data is correct and sufficient. If this is not the case, then OpenPhase will show 
an error message and expect the user to correct the input file independently. 


9.4.1 Step by Step Through the ProjectInput File 


This section will walk through a default ProjectInput file line by line, explaining 
exactly what each line does and what it represents. 


RunTimeControl 


The input file begins with the RunTimeControl module input. This controls the 
simulation time, the frequency of outputs, and whether the user wants to start a new 
simulation or restart an existing simulation from a specific time step. Figure 9.1 
shows all the parameters related to the Run TimeControl module. The following 
gives a description of each parameter: 


— SimTtl 
The user can specify a title for the simulation. 


Fig. 9.1 The RunTimeControl section in the ProjectInput file 


9.4 The Input File Structure 


— nSteps 
The number of time steps for the entire simulation. To calculate the actual 
time, you have to multiply this number by At. 


— FTime, STime 
Sets parameters for the frequency of data output to the disk (FT ime) and to 
the screen (STime). 


— LUnits, TUnits, MUnits 
All the values are in the SI system of units, which uses the meter, second, and 
kilogram as base units. 


— dt 
This is the initial time step (time discretization) At. Note that some modules 
might use an internal time-stepping scheme. 


— nOMP 

The number of OpenMP threads can be specified if you want to run your 
simulation in a multithreading mode. The number 1 indicates no active 
multithreading. Make sure to understand how OpenMP works before running 
your simulation using multithreading because if you increase it unwisely, it 
can slow your simulation instead of making it faster. 


— Restrt 
OpenPhase gives the user the option of either starting a new simulation (No) 
or restarting from the existing previously saved output (Yes). 


— tStart 
Suppose you indicate (Yes) in the previous option. In that case, you have to 
specify the time step from which you want to restart your simulation. 
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— tRstrt 
Set parameter for the frequency of raw-data output to the disk for restarting. 


Settings 


The Settings input controls the size and numerical resolution of the simulation 
domain, and the user must specify the grid spacing and the number of grid points 
along each dimension. The grid spacing is one of the most critical numerical 
parameters because it defines the resolution of the simulation. In practice, high 
resolution means more time, so it is usually necessary to find a middle ground 
between grid resolution and calculation time. Other information linked to the phases 
is contained inside the Settings parameters, including the thermodynamic phase 
of each phase field and various additional information about the phases, such as 
their aggregate state (solid, liquid, or gas) and the chemical components, among 
other things, as shown in Fig. 9.2. 


— Nx, Ny, Nz 
The simulation system size in the x, y, and z directions given in grid points. 


Fig. 9.2 The Settings section in the ProjectInput file 
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— dx 
Sets the grid spacing Ax. It is important to note that this value is the same for 
all three spatial dimensions. 


— iWidth 
The interface width ņ of the phase field given in grid points. 


— Phase 
The name of the phase; the running index refers to the thermodynamic phase 
for a given phase field. 


— State 
The aggregate state of the phase (solid, liquid, or gas). 


— Comp 
The chemical component name. 


— Nvariants 
The number of crystallographic (symmetry/translation/...) variants for a given 
phase. 


BoundaryConditions 


The boundary conditions of the simulation domain have to be set individually for 
every simulations domain side in the Project Input file, as shown in Fig. 9.3. 
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Fig. 9.3 The BoundaryConditions section in the ProjectInput file 


— BCNO0X, BCNX, BCNOY, BCNY, BCN0Z, BCNZ 
These keywords define the six individual sides of the calculation domain. 


In OpenPhase The following boundary conditions are available: 


* Periodic: Periodic boundary conditions. 

* NoFlux: Adiabatic boundary conditions, resulting in zero first-order derivative at 
the boundary. 

* Free: Boundary conditions resembling a free surface, resulting in a continuous 
gradient across the boundary. 

* Fixed: Fixed (Dirichlet or Neumann) boundary conditions. The value should be 
set by the user. 
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Chapter 10 A) 
Tutorial 2: OpenPhase Examples TCA 


10.1 Normal Grain Growth 


Grain growth is a microstructure-evolution phenomenon that occurs in polycrys- 
talline materials after solidification. In this process, some grains are growing at 
the cost of the shrinkage of other grains or crystallites of a specific orientation. 
Consequently, the average grain size increases, and the number of grains in the 
material decreases. 

Grain growth is curvature-driven grain-boundary motion between different 
crystallites with different orientations. The energy of these grain boundaries adds 
excess free energy to the system. Therefore, the reduction of the grain-boundary 
energy is the driving force for the movement of the grain boundaries, and the system 
evolves towards minimizing the total free energy. 

As a result, the material's microstructure proceeds to reduce the total interface 
area. The kinetics of grain boundaries influence the interfacial migration; the energy 
and mobility of these boundaries control how the grain structure evolves. 


10.1.1 Simulation Example 


The OpenPhase distribution comes with the NormalGG example, which demon- 
strates the simulation of grain growth. The driving force is curvature, which is 
the only contribution to the microstructure evolution. Other contributions, such as 
temperature and concentration, can easily be coupled. 

The example shows normal grain growth with isotropic interface energy and 
mobility, and the user can modify the example to account for the anisotropy of both 
contributions. 
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Table 10.1 Simulation 


Parameter Value Units 
pare Grid size 1013 Cells 
Grid spacing Ax Ix 1076 (m 
Time step At 1x10 |$ 
Interface width n 5Ax m 
Interface energy o 0.24 zs 


Interface mobility u | 1 x 1077 


The microstructure is generated using Voronoi tessellation. The initial num- 
ber of grains can be changed directly by the user in the example file. The 
ProjectInput.opi file contains all the input parameters related to the sim- 
ulation, and the user can modify it for different applications. In this simulation, 
random crystallographic orientations are assigned to all grains. Periodic boundary 
conditions are introduced in all directions. The other simulation parameters are 
listed in Table 10.1. 

The reader can find the source code of this example in the OpenPhase distribution 
in the examples folder. The source code is also included in this chapter (Listing 1). 
Warning: the code presented here is only for illustration and is subject to changes in 
the future. The reader is therefore referred to the OpenPhase website for the actual 
documentation. 


Step by Step Through the Source Code 


Lines 1-23: Importing different classes from the OpenPhase library that are used 
in the simulation of this example. 

Lines 26-28: This function sets the initial microstructure as shown in Fig. 10.1b. 
In this illustrative example, 200 grains have been initialized with the same 
thermodynamic phase but different phase-field indices based on the Voronoi 
tessellation technique implemented in the Initializations class. The 
grains have a random crystallographic orientation, indicated by the RGB colors 
in Fig. 10.1a. 

Line 32: This function sets the interface properties: mobility and energy. In this 
example, isotropic interface energy and mobility have been used, as shown in 
Fig. 10.2. However, the user can use different anisotropy model forces for the 
interface energy and mobility in the ProjectInput file. 

Line 33: This function calculates the interface-curvature-related driving force. It 
includes the multi-junction energy as indicated in Eq. (6.1). 

Line 34: This function limits the phase-field increments for all present phase- 
field pairs so that the actual phase-field values are within their natural limits of 0 
and 1. 

Line 35: This function merges the increments into the phase fields. It finalizes 
the phase-field calculations by setting the boundary conditions, calculating the 
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Normal Grain Growth 


#include "Settings.h" 

#include "RunTimeControl.h" 
#include "DoubleObstacle.h" 
#include "PhaseField.h" 

#include "Initializations.h" 
#include "BoundaryConditions.h" 
#include "InterfaceProperties.h" 
#include "Tools/TimeInfo.h" 
#include "Info.h" 


using namespace std; 
using namespace openphase; 


int main() 
{ 
string InputFile = "ProjectInput.opi"; 
Settings OPSettings; 
RunTimeControl RTC (OPSettings) ; 
PhaseField Phi (OPSettings) ; 
DoubleObstacle DO (OPSettings) ; 
InterfaceProperties IP(OPSettings); 
BoundaryConditions BC(OPSettings); 
TimeInfo Timer(OPSettings, "Execution Time 


> Statistics"); 


//generating initial grain structure using Voronoi algorithm 
int number of grains - 200; 

size t GrainsPhase - 0; 
Initializations::VoronoiTesselation(Phi, BC, OPSettings, 

— number of grains, GrainsPhase); 


for(RTC-tStep = RTC-tStart; RVC tStep <= RIC nSteps; 
> RTC.IncrementTimeStep () ) 


{ 
TESSUTI) 
DO.CalculatePhaseFieldIncrements (Phi, IP); 
Phi.NormalizeIncrements (BC, RTC.dt); 
Phi.MergeIncrements (BC, RTC.dt); 
if (RTC.WriteVTK() 
1 
Phi.WriteVTK(RTC.tStep, OPSettings); 
} 
if (RTC.WriteToScreen () 
{ 
double I En = DO.AverageEnergyDensity (Phi, IP); 
std::string message = Info::GetStandard("Interface energy 
< density", to string(I En)); 
Info::WriteTimeStep(RTC, message); 
Timer.PrintWallClockSummary(); 
} 
} 
return 0; 


Listing 1: Source code for the NormalGG example in OpenPhase 
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© 
(a) (b) 


Fig. 10.1 Initial microstructure, displayed in different modes: (a) RGB colors; (b) as phase-field 
index, where the diffuse interface region is displayed in light grey 


$MobilityModel Interface mobility model : IS © 


Fig. 10.2 The InterfaceProperties module in the ProjectInput file 


derivatives, setting the flags that mark interfaces, collecting the volume for each 
phase field, and calculating the thermodynamic phase fractions. 

Line 39: This function writes an output file in VTK format, which contains the 
information of the phase field (such as the interfaces, the junctions, and the phase 
fractions of the different phases). III ParaView can be used to visualize the VTK 
files, and the output results can be seen in Fig. 10.3. 

Lines 41-47: This code block prints different simulation information on the 
screen. 


Results 


Figure 10.3 shows the microstructural evolution for the given parameters, in which 
the grain boundaries appear as grey regions separating individual grains. 
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(b) 


Fig. 10.3 Evolution of the microstructure at different times f (in time steps). (a) t = 500. (b) 
t = 1500. (c) t = 3000 


10.2 Dendritic Solidification 


Phase-field simulation can be used to model the solidification of metallic alloys. 
Interface velocities are well represented by a Gibbs-Thomson relation. This is 
known as the Stefan problem when linked to heat and solute redistribution at the 
interface and long-range transport in bulk phases. Interfacial energy anisotropy is 
weak and can be considered a perturbation of a spherical Wulff shape. Even though 
the solid and liquid have different densities, no significant stresses are developed 
at the interface during growth. In comparison to solute diffusion in the melt, heat 
conduction in both phases is almost the same, and solute diffusion in the solid can be 
ignored [1]. In general, the evolution of the dendrite interface area is characterized 
by a rapid increase during the growth of side branches, and this is followed by a 
decrease due to capillarity and coalescence of adjacent structures. 


10.2.1 Simulation Example 


The OpenPhase distribution provides the So1idificationFeC example, which 
is a simulation of dendritic solidification.! The code presented here is written 
considering a 2D case. However, it can be modified to account for the 3D case. In 
this example, two different simulations have been considered: 2D equiaxed dendritic 
solidification for a single grain and for multiple grains with different orientations 
growing from the melt. 

The reader can find the source code for this example in the OpenPhase dis- 
tribution in the examples folder. The source code is also included in the chapter 


' The code presented here was developed in collaboration by a Chinese exchange student, MSc. 
Jianghai Cao from Chongqing University, China, during his interim at RUB. 
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(Listing 2). Warning: the code presented here is only for illustration and is subject 


to changes in the future. Therefore, the reader 
for the actual documentation. 


#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
using namespace std; 

using namespace openphase; 
fee <<< The Main 


int main(int argc, char «argv[l) 


{ 


"Settings.h" 
"RunTimeControl.h" 
"InterfaceProperties.h" 
"DoubleObstacle.h" 
"PhaseField.h" 
"DrivingForce.h" 
"Composition. h" 
"Temperature. h" 


"BoundaryConditions.h" 
"Pnitiali zations. n" 
"Tools/TimeInfo.h" 


>35 


Settings 
OPSettings.ReadInput(); 
RunTimeControl 
PhaseField 
DoubleObstacle 
nterfaceProperties 
EquilibriumPartitionDiffusionBinary 
Composition 
Temperature 
DrivingForce 
BoundaryConditions 
TimeInfo 


& Statistics"); 


idx0 
i.FieldsStatistics[idx0].State 
idx1 Phi.PlantGrainNucleus(1, 


OPSettings.Nz/2); 


Cx.SetInitial(BC, Phi); 
Tx.SetInitial (BC); 
DF.SetDiffusionCoefficients (Phi, 


The Time Loop 


—  RTC.IncrementTimeStep () ) 


{ 


"Set (Phi); 


p 


1j 


.GetDrivingForce (Phi, 
.Average(Phi, BC); 


(CoS WESS 


hi.NormalizeIncrements (BC, RTC. 
F.Solve(Phi, Cx, Tx, 
LEESE, Phi, RTC. dt); 

hi.MergeIncrements (BC, RTC.dt); 


= 
4 


p 


Initializations::Single(Phi, 
AggregateStates::Liquid; 


Tx); 


.CalculateInterfaceMobility (Phi, 
.CalculatePhaseFieldIncrements (Phi, 


.MergePhaseFieldIncrements (Phi, 


is referred to the OpenPhase website 


"EquilibriumPartitionDiffusionBinary.h" 


"Tools/MicrostructureAnalysis.h" 


ook eoe ee e ee / 


OPSettings; 


RTC(OPSettings); 
Phi (OPSettings); 
DO(OPSettings); 
P(OPSettings 


( DE 
DF(OPSettings); 
Cx(OPSettings); 
Tx(OPSettings) 
dG (OPSettings); 
BC(OPSettings); 


Timer (OPSettings, 


i 


"Execution Time 


0, BC, OPSettings) ; 


OPSettings.Nx/2, OPSettings.Ny/2, 


i.FieldsStatistics[idxl].State = AggregateStates: :Solid; 


i 


RUNG CSEart; RTE- CSCEp <= KING nsuempsm 


(che. 184, BG; 
IP); 


TP); 


dG) ; 


TES 
ati; 


BC, RTC.dt); 
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56 if (RTC.WriteVTK() ) 


57 { 

58 Write data in VTK format 

59 Phi.WriteVTK(RTC.tStep,OPSettings); 

60 Cx.WriteVTK(RTC.tStep,OPSettings); 

61 Tx.WriteVTK(RTC.tStep,OPSettings); 

62 IP.WriteVTK(Phi,RTC.tStep); 

63 } 

64 Output to screen 

65 if (RTC.WriteToScreen()) 

66 { 

67 double I En - DO.AverageEnergyDensity(Phi, IP); 

68 std::string message = Info::GetStandard("Interface energy 
< density", to string(I En)); 

69 Info::WriteTimeStep(RTC, message); 

70 dG.PrintDiagnostics(); 

71 Phi.PrintVolumeFractions(); 

72 Timer.PrintWallClockSummary () ; 

7 n 

74 } end time loop 

75 return 0; 


Listing 2: Source code for the SolidificationFeC example in OpenPhase 


Step by Step Through the Source Code 

1: Equiaxed Dendritic Solidification for a Single Grain 

Lines 1-30: Importing different classes from the OpenPhase library that are used 
in the simulation of this example. 

Lines 32-35: These functions set the initial microstructure as illustrated in 
Fig. 10.4a. There is a matrix with a single phase representing the liquid, which 
has a phase index equal to 0, and we plant a grain nucleus in the middle of our 
simulation domain, which represents the solid phase and has a phase index equal 
to 1. 

Lines 37-38: These functions set the initial composition, and the temperature 
and diffusion coefficients for the system. They take their inputs from the 
ProjectInput file. 

Line 44: This function sets the interface properties: mobility and energy. In this 
example, anisotropic driving forces for the interface energy and mobility have 
been set in the ProjectInput file. 

Line 45: This function calculates concentration-dependent mobility. 

Line 46: This function calculates the interface-curvature-related driving force. It 
includes the triple-junction energy as indicated in Eq. (6.1). 

Line 47: This function calculates the driving force for each point. 

Line 48: This function averages the driving forces across the interface. 
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Line 50: This function calculates a time-independent phase-field increment by 
converting the local driving forces. An additional noise term can be applied, as 
can a stability mechanism to strengthen the numerical phase-field profile. 

Line 51: This function limits phase-field increments for all present phase-field 
pairs so that the actual phase-field values are within their natural limits of 0 and 
1. 

Line 52: This function calculates the change of total concentrations in one time 
step taking into account cross terms. 

Line 53: This function sets the temperature according to the cooling rate; it 
applies equal temperature increments to all points. 

Line 54: This function merges the increments into the phase fields. It finalizes 
the phase-field calculations by setting the boundary conditions, calculating the 
derivatives, setting the flags that mark interfaces, collecting the volume of each 
phase field, and calculating the thermodynamic phase fractions. 

Lines 57-65: This code block writes different output files in VTK format, which 
contains the information of the phase field (such as the interfaces, the junctions, 
and the phase fractions of the different phases), composition, temperature, and 
interface properties. ZII ParaView can be used to visualize the VTK files. 

Lines 66—74: This code block prints different simulation information on the 
screen. 


2: Equiaxed Dendritic Solidification for Multiple Grains with Different Orien- 

tations 

For this case, the code stays the same as discussed above with a little modification 

to the initialization part, in which line 34 can be replaced by the following line: 
This function plants random nucleation sites in the simulation domain with 

different orientations. The number of particles can be chosen by replacing 

NoOfParticles with an integer number. 


Results 


Figures 10.4 and 10.5 show that side branches emerge as the dendritic tips extend 
into the undercooled melt, resulting in a complicated solid-liquid interface morphol- 
ogy. Once you reach a certain distance behind the dendritic tip, a slower evolution 
begins to take place, which occurs at or around the point of phase equilibrium. 
Capillarity and delayed solidification play a significant role in determining the 
interface dynamics at this stage. 
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Fig. 10.4 Equiaxed dendritic solidification in 2D. (a) t = 0.05s. (b) t = 0.25s. (c) t = 0.75s. 
(d): = 1.5s 
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Fig. 10.5 Equiaxed dendritic solidification of multiple grains with different orientations. (a) t = 
0.05s. (b) t = 0.25s. (c) t = 0.75s. (d) t = 1.5s 
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Appendix A 


A. Simple Dendrite Code 


This section presents the “100-line C++ code" that we use to teach the students 
how to simulate a dendritic growth similar to the example illustrated in Chap. 1. A 
good exercise for the reader is to choose suitable parameters and correctly repeat 
the simulation given in “Example—Dendritic Growth" in Chap. 1. 


<cmath> 
<fstream> 
<sstream> 
<iostream> 
«cstdlib» 
«complex» 
«random» 


1 #include 
2 #include 
3 #include 
4 #include 
5 #include 
6 #include 
7 #include 
8 

9 


using namespace std; 


11 // Discretization parameters 


12 const int Nx - 300; 

13 const int Ny = 300; 

14 const int Nt - 

15 const int BCells = 1; 

16 double dx ="/0503% 
17 double dy = 0.03; 
18 double dt = 2.0e-4; 


20 // Kobayashi's parameters 


21 double epsilon = 0.010; 
22 double tau = 3.0e-4; 
23 double alpha = 0.7} 
24 double Gamma = 10.0; 
25 double delta = 0.04; 
26 double theta0 - M PI/2; 
27 int sym z 8; 
28 double K = _ 2.09 

> (no-dimension) 
29 double TO - 0.0; 
30 double Te = q.0; 

—  (no-dimension) 
31 int RandomSeed = 125; 
32 double ampl = 0.02; 
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// 
// 
Af 
// 


fi. 
rd d 


"74 


/ Anisotropy in 
^ Rotation of the crystal orientation 


// Discretization x-direction 
// Discretization y-direction 
1000000000;// Number of time steps 


Number of boundary cells 
Grid spacing in x-direction 
Grid spacing in y-direction 
Size of time step [s] 


[m] 
[m] 


(not exactly his..) 


Gradient energy coefficient 
Inverse of interface mobility [s] 
Coefficient of driving force 
Coefficient of driving force 
(0,1) 


Number of symmetry fold 
Referrers to latent heat 


Initial temperature 
Equilibrium temperature 


Defines the set of random numbers 
Amplitude of noise 
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// Misc parameters 


int tout = 100; // Output distance in time steps 

double t = 0.0; // Simulation start time [s] 

double Radius = 0.3; // Initial radius of spherical grain 
double PhiPrecision = 1.e-9; // Phase-field cut off 

double TempPrecision = 1.e-3; // Precision of Gauss-Seidel algorithm 
int TempMaxIter - 10000; // Maximum number of Gauss-Seidel 


— iterations 


// Phase-field and Temperature storage 


const int Mx = Nx + 2 * BCells; // Shorthand for entire storage loop 
const int My = Ny + 2 + BCells; // Shorthand for entire storage loop 
double Phi [Mx] [My] ; // Phase-field storage 

double PhiDot [Mx] [My]; // dPhi dt storage 

double Temp [Mx] [My] ; // Temperature storage 

double Eps [Mx] [My] ; // Laplace parameter 

double EpsP [Mx] [My] ; // Derivative of Laplace parameter 
double Noise [Mx] [My] ; // Storage for noise 


/** Calculates Inclination angel 


* \param d position on grid in x-direction 
* \param j position on grid in y-direction 
* \return angel theta in [0,pi] 
*x/ 
double CalculateTheta(int i, int j) 
{ 
double dPhix = (Phi[i+1][j ] - Phi[i-1] [j 1)/(2«8x); 
double dPhiY = (Phi[i ]l[j«1] - Phili ] [j-1])/(2x«dy); 


double dPhiNorm = sqrt( pow(dPhiX, 2) + pow(dPhiY, 2)); 


double theta - 0.0; 

if (dPhiNorm » 1.e-10) 

{ 
double NormX = dPhiX/dPhiNorm; 
double NormY  - dPhiY/dPhiNorm; 
theta += atan2(NormY, NormX); 

} 


return theta; 


li 


/** Anisotropy function 


* \param theta Angel between x-axis and interface normal 
* \return Amplitude of Laplace term 
*x/ 


double sigma(double theta) 


{ 


double value = 1 + delta * cos( sym « (theta - theta0)); 
return value; 


)3 


/*»« First derivative of anisotropy function 


* \param theta Angel between x-axis and interface normal 
* \return Amplitude of Laplace term 
*x/ 


double sigmaP (double theta) 


{ 


double value = - delta * sym « sin( sym » (theta - theta0)); 
return value; 


li 


/** Initializes a grain of a certain radius in a supercooled melt 
* \param Radius Radius of grain 


Appendix A 

99 *x/ 
100 void InitializeSupercooledGrain (double Radius) 
w1 ( 
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// Initialization 

const double x0 = Nx/2 + dx; 

const double y0 = 0; 

for (int i = 0; i < Nx + 2 « BCells; i++) 
for (int j 0; j < Ny + 2 x BCells; j++) 


{ 


// Initialize the phase Field 


double r = sqrt (pow(ixdx-x0,2) + pow(j«dy-y0,2)); 


if ( r < Radius) 

{ 
Phi [i] [j] = 
PhiDot [i] [j] = 


or 
oS) 


} 


else 


{ 

Phi [i] [j] = 0.0; 
PhiDot [i] [j] = 0.0; 
} 
// Initialize temperature 


Temp [i] [3] = TO; 


} 


/x* Calculate Laplacian of the phase-field 
* \param i Position on grid in x-direction 
* \param j Position on grid in y-direction 
* \return Laplacian of phase-field 
*x/ 

double LaplacePhi(int i, int j) 


{ 


double Laplacian = 0.0; 
( 


Laplacian += 2./3.«(Phi[i«1][j ] - 2.0«Phililljl 


=> 1)/(dx+dx); 


Laplacian += 2./3.«(Phi[i ][j+1] - 2.0«Philil [j] 


=> ]llij-i1D/(dy«dy); 


Laplacian += 1./3.«(Phi[i«1][j«1] + Phi[i+1] [j-1] 
>  Phi[i-1][j-1] - 2«Phili] [j1)/ (pow (sqrt (2.)*dx,2)); 


return Laplacian; 


} 


/** Calculate Laplacian of the temperature 
* \param i Position on grid in x-direction 
* \param j Position on grid in y-direction 
* \return Laplacian of temperature 
**/ 

double LaplaceTemp(int i, int j) 


{ 


double Laplacian = 0.0; 


Laplacian += 2./3.«(Temp[i«1][j ] - 2.0«Temp[il [j] 


> ])/(dx«dx); 


Laplacian += 2./3.«(Temp[i ][j+1] - 2.0«Temp[i] [j] 


> ]llij-11)/(dy«dy) 


+ Phili-i][j 


+ Phi[i-1] [j+1] 


+ Temp [i-1] [j 
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Laplacian += 1./3.*(Temp[i+1] [j+1] + Temp[i+1] [j-1] + Temp[i-1] [j+1] + 


> Temp[i-1] [j-1] - 4«Temp[i] [j1) / (pow (sqrt (2.)*dx,2)); 


return Laplacian; 


} 


/** Writes output to file 
* \param tStep time step 
**/ 

void WriteToFile(int tStep) 


{ 
{ 


stringstream outbufer; 


outbufer -- 


"# vtk DataFile Version 3.0\n"; 


outbufer << "PhaseField\n"; 
outbufer << "ASCII\n"; 
outbufer << "DATASET STRUCTURED _GRID\n"; 
outbufer << "DIMENSIONS " << Nx << " " << Ny << " " 
outbufer << "POINTS " << Nx«Ny«l1 << " int\n"; 
for(int j = 1; j < Nx+1; ++j) 
for(int i = 1; i < Ny+l; ++i) 

outbufer sad ee "o" åa j xc feg ee "Ap 
outbufer << "Mn"; 
outbufer << "POINT DATA " << Nxx+Nyxl << "Xm"; 
outbufer << "SCALARS PhaseFields double 1" << " Mn"; 
outbufer << "LOOKUP TABLE default" << " \n"; 
for(int j = 1; j < Nx+1; ++j) 
for(int i = 1; i « Ny+1; ++i) 

outbufer << double(Phi[i][j]) << " n"; 
outbufer << "SCALARS Temperature double 1" << " \n"; 
outbufer << "LOOKUP TABLE default" << " Wn"; 
for(int j = 1; j < Nx+1; ++j) 
for(int i = 1; i « Ny+1; ++i) 

outbufer << double (Temp [i] [j]) << " \n"; 


} 


stringstream filename; 


filename << "PhaseField_" << tStep << ".vtk"; 


string FileName = filename.str(); 
ofstream vtk file(FileName.c str()); 
vtk file -- outbufer.rdbuf(); 

vtk file.close(); 


/** Main Kobayashi code +*/ 
int main() 


{ 


// initialize random seed 

srand (RandomSeed) ; 
InitializeSupercooledGrain (Radius) ; 
WriteToFile(0); 


// Main time loop 


for 


{ 


(int tStep = 1; tStep < Nt; tStep++ 


// Calculate gradient coefficient 


for (int i = BCells; i < Nx + BCells; i++) 
for (int j = BCells; j < Ny + BCells; j++) 
{ 
if ((Phi[i] [j] > PhiPrecision) and (Phi[i] [j] < 
— PhiPrecision) ) 
{ 
double theta = CalculateTheta(i,j); 
Eps [i] [j] = epsilon « sigma (theta); 
EpsP[i] [j] = epsilon +» sigmaP (theta); 
} 
else 
{ 
Eps [i] [j] = epsilon; 
EpsP[i] [j] = 0.0; 


<< I xe 


Ts 
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"An"; 
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226 
227 
228 
229 
230 
231 
232 
233 
234 
235 
236 
237 
238 


239 
240 


241 


242 
243 
244 
245 


246 


247 


248 


249 
250 
251 
252 
253 
254 
255 
256 
257 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 
278 
279 
280 
281 
282 
283 
284 
285 


} 

double noise = Phi[i] [j] * (1 -Philil[jl) * ampl; 
noise «= Noise[i] [j]; 

Eps [i] [j] += noise; 


) 


// Calculate PhiDot 
for (int i = BCells; i < Nx + BCells; i++) 
for (int j = BCells; j < Ny + BCells; j++) 


{ 
PhiDot [i] [j] = 0.0; 
// Calculate driving force m 
if ((Phi[i][j] > PhiPrecision) and (Phi[i][j] < 1. - 
> PhiPrecision) ) 
{ 
double m = (alpha/M PI) +» atan(Gammax (Te - 
>  Temp[ilIj1)); 
PhiDot [i] [j] += Phili][jl«(1.0 - Phil[il[jl)«(Philil[j] - 
— 0.5 « m); 
) 
// Calculate Laplacian-term 
PhiDot [i] [j] -= (Eps[i+1] [j] * EpsP[i+1] [j] » (Phi[i+1] [j+1] 
<> Phi[i+1] [j-1]) 
- Eps[i-1] [j] +» EpsP[i-1][j] * (Phi[i-1] [j+1] 
>  Phil[i-1][j-11))/ (4«dx-«dy); 
PhiDot [i] [j] += (Eps[i] [j+1] * EpsP[i] [j+1] « (Phi[i+1] [j+1] 
— Phi[i-1] [j+1]) 
- Eps[i] [j-1] * EpsP[i] [j-1] * (Phi[i+1] [j-1] 
— Phi[i-1] Pn D 
double dPhix = (Phi[i+1] [j] - Phi[i-1 )/(2+dx) ; 
double dPhiY = (Philil[j«1] - Phil[i] m Dos 
double dEpsX = (Eps[i+1] [j - Eps [i-1 1)/(2«àx) ; 
double dEpsY = (Eps[il[j«1] - Eps[i] 5-11) /( 2+dy) ; 
PhiDot [i] [j] += 2«Eps[i][j]* (dPhiX«dEpsX + dJPhiYsdEpsY) ; 
PhiDot [i] [j] += pow (Eps [i] [j],2) » LaplacePhi(i,j); 
PhiDot [i] [j] /= tau; 
} 


// Update Phase-field 
for (int i = BCells; i < Nx + BCells; i++) 
for (int j BCells; j « Ny + BCells; j++) 


{ 
if (PhiDot [i] [j] != 0.0) 
{ 
// Update phase field 
Phi [i] [j] += dt * PhiDot [i] [j]; 
// Limit phase field 
if (Phi[i] [j] < PhiPrecision) Phi[i] [j] = 0.0; 
else if (Phi[i] [j] > 1 - PhiPrecision) Phi[i] [j] = 1.0; 
} 
} 


// Apply adiabatic boundary conditions 
for (int i = 0; i < Nx + 2«BCells; i++) 


{ 
Phi[i] [0 ] = Phi[i] [1 lis. 
Phi [i] [Ny] = Phili] [Ny-1]; 
} 
for (int j = 0; j < Ny + 2*BCells; j++) 
{ 
Phi[0 ] [j] = Phil[1 1s 
Phi [Nx] [j] = Phi [Nx-1] [j]; 
} 


// Calculate TempDot (implicit scheme with Gauss-Seidel) 
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double MaxDTemp - 
int iterations - 
do 
{ 
MaxDTemp = 0.0; 
// Save current Temp 
for (int i = BCells; 
for (int j = BCells; 
{ 
double TempOld = 
double aij - (- 
double bij - 
double neightbou 
neightbours += d 
neightbours += d 
neightbours «- d 
neightbours += d 


Temp [i] [j] = 1./ 
MaxDTemp = max (M 


iterations++; 


} 
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// Maximum of Temp change in on iteration 


// Iteration counter 


erature in TOld 


// Apply adiabatic boundary conditions 


for (int i = 0; 2. < 
{ 
Temp [i] [0 ] 
Temp [i] [Ny 


) 


for (int j = 0; j < 
( 
Temp [0 JLH -= 
Temp[Nx ][j] = 


} 
} 


while 


((MaxDTemp > TempP 


// Calculate calculation 
t += dt; //increase the 
// Write into VTK File 
if( tStep%tOut == 0) 
{ 
WriteToFile(tStep) ; 
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